show code
knitr::opts_chunk$set(echo = TRUE)Nucleation Experiment
After first submission
Data and R codes available: github, then it will be submited to zenodo
Inclusion of species diversity of recruiting plants: item I.C, with rarefaction and Proportion of Biotic/Abiotic Seed Dispersal and Planted/Unplanted species.
Table with statistical details from the models (file
"itatinga_NUC_ModelSummaries.xlsx")Fixed LiDAR figure (SE of the mean)
Fixed plots with empty recruits
Considering only small recruits in species richness rarefaction curves
Analyses separating Baccharis + Vernonanthura
Setup environment
1) R Environment
show code
path <- getwd()
version _
platform x86_64-w64-mingw32
arch x86_64
os mingw32
crt ucrt
system x86_64, mingw32
status
major 4
minor 3.3
year 2024
month 02
day 29
svn rev 86002
language R
version.string R version 4.3.3 (2024-02-29 ucrt)
nickname Angel Food Cake show code
if(!require(pacman)){install.packages("pacman")}Carregando pacotes exigidos: pacmanshow code
pkg_list <-c("tidyverse", "lme4", "DHARMa", "emmeans", "multcomp", "iNEXT",
"car", "performance", "cowplot", "MuMIn", "ggpattern")
pacman::p_load(pkg_list, character.only = TRUE)2) Read data
trees_clean: data of trees (planted/not-planted)
show code
trees_clean <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_trees.xlsx"))
trees_clean <- trees_clean |>
dplyr::mutate(scientificName = dplyr::if_else(is.na(scientificName), "Unclassified", scientificName))grass_cover: grass cover per plot (%)
show code
grass_cover <- readxl::read_xlsx(paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))lidar: Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).
show code
lidar2022_wide <- readxl::read_xlsx(
paste0(path, "/tab_metrics_LiDAR_2022_v2.xlsx")) |>
dplyr::select(Trat:LAIunder_SD) |>
tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
ano = 2022) |>
dplyr::select(-Perc_cobertura) New names:
• `` -> `...1`Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].show code
lidar2024_wide <- readxl::read_xlsx(
paste0(path, "/tab_metrics_LiDAR_2024_v2.xlsx"))|>
dplyr::select(Trat:LAIunder_SD) |>
tidyr::separate(Trat, into = c("Bloco", "Perc_cobertura","Tratamento"), sep = "_") |>
dplyr::mutate(Tratamento = dplyr::if_else(is.na(Tratamento), "controle", Tratamento),
Tratamento = paste(Perc_cobertura, Tratamento, sep="_"),
ano = 2024) |>
dplyr::select(-Perc_cobertura)New names:
• `` -> `...1`Warning: Expected 3 pieces. Missing pieces filled with `NA` in 5 rows [5, 10, 15, 20,
25].show code
lidar_wide <- dplyr::bind_rows(lidar2022_wide, lidar2024_wide) |> # 40 x 12
dplyr::mutate(
Tratamento = factor(Tratamento,
levels = c("25_tabuleiro", "25_faixa", "50_tabuleiro", "50_faixa", "100_controle"),
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")))Table - List of sampled species (planted and recruits)
two tables:
show code
trees_clean |> dplyr::select(RP, family, scientificName) |>
unique()|>
dplyr::mutate(valor=TRUE)|>
tidyr::pivot_wider(names_from = RP, values_from = valor, values_fill = FALSE) |>
dplyr::arrange(desc(planted), family, scientificName) # |> writexl::write_xlsx("tab1.xlsx")# A tibble: 108 × 4
family scientificName planted regenerant
<chr> <chr> <lgl> <lgl>
1 Anacardiaceae Astronium graveolens TRUE FALSE
2 Annonaceae Annona cacans TRUE TRUE
3 Bignoniaceae Handroanthus impetiginosus TRUE FALSE
4 Bignoniaceae Handroanthus ochraceus TRUE TRUE
5 Boraginaceae Cordia sellowiana TRUE TRUE
6 Euphorbiaceae Croton floribundus TRUE TRUE
7 Euphorbiaceae Croton urucurana TRUE FALSE
8 Euphorbiaceae Mabea fistulifera TRUE TRUE
9 Fabaceae Anadenanthera colubrina var. cebil TRUE TRUE
10 Fabaceae Apuleia leiocarpa TRUE FALSE
# ℹ 98 more rowsshow code
trees_clean |> dplyr::select(RP, family, scientificName) |>
unique() |> with(table(RP))RP
planted regenerant
37 86 I. Results - Field data
A) Survival of planted trees
(i) Organizing data
show code
trees_clean |> with(table(RP, scientificName)) |>
as.data.frame() |>
dplyr::filter(Freq>0) |> # writexl::write_xlsx("tab_spp_plantio.xlsx")
dplyr::filter(RP=="planted") RP scientificName Freq
1 planted Anadenanthera colubrina var. cebil 204
2 planted Annona cacans 4
3 planted Apuleia leiocarpa 11
4 planted Astronium graveolens 5
5 planted Cariniana estrellensis 2
6 planted Cecropia pachystachya 117
7 planted Cedrela fissilis 28
8 planted Ceiba speciosa 54
9 planted Citharexylum myrianthum 16
10 planted Copaifera langsdorffii 24
11 planted Cordia sellowiana 10
12 planted Croton floribundus 831
13 planted Croton urucurana 726
14 planted Enterolobium contortisiliquum 77
15 planted Erythrina mulungu 32
16 planted Eugenia uniflora 9
17 planted Ficus guaranitica 82
18 planted Gallesia integrifolia 50
19 planted Genipa americana 3
20 planted Guazuma ulmifolia 323
21 planted Handroanthus impetiginosus 30
22 planted Handroanthus ochraceus 20
23 planted Hymenaea courbaril 12
24 planted Inga edulis 85
25 planted Inga vera 741
26 planted Lafoensia pacari 34
27 planted Luehea divaricata 48
28 planted Mabea fistulifera 20
29 planted Machaerium acutifolium 17
30 planted Myroxylon peruiferum 2
31 planted Peltophorum dubium 29
32 planted Platypodium elegans 7
33 planted Psidium myrtoides 2
34 planted Pterogyne nitens 3
35 planted Senegalia polyphylla 474
36 planted Senna macranthera 1
37 planted Unclassified 372Reading scientific names
show code
survival_planted2024 <- trees_clean |> # 5416 stems
dplyr::filter(ano==2024)|>
dplyr::select(bloco:arv, Hm, scientificName, RP) |> # 3,199 stems
unique() |> # 1487
dplyr::filter(RP=="planted") |> # 1135
dplyr::group_by(tratamento,bloco, parcela)|>
dplyr::summarise(n = dplyr::n(), .groups = "drop")|>
dplyr::arrange(desc(n))(ii) Survival analysis
24 seedlings per plot
show code
survival_planted2024 |>
with(table(bloco, parcela, by=tratamento)), , by = 100
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = F25
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = F50
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = N25
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1
, , by = N50
parcela
bloco 1 2 3 4
1 1 1 1 1
2 1 1 1 1
3 1 1 1 1
4 1 1 1 1
5 1 1 1 1show code
survival_planted2024 |> dplyr::arrange(tratamento, bloco, parcela)# A tibble: 100 × 4
tratamento bloco parcela n
<chr> <chr> <chr> <int>
1 100 1 1 17
2 100 1 2 15
3 100 1 3 12
4 100 1 4 19
5 100 2 1 11
6 100 2 2 11
7 100 2 3 17
8 100 2 4 13
9 100 3 1 13
10 100 3 2 16
# ℹ 90 more rowsshow code
mod_surv1 <- lme4::glmer(cbind(n, 24 - n) ~ tratamento + (1|bloco), data = survival_planted2024, family = binomial(link = "logit"))
mod_surv_null <- lme4::glmer(cbind(n, 24 - n) ~ 0 +(1|bloco), data=survival_planted2024, family = binomial(link = "logit"))
AIC(mod_surv1, mod_surv_null) df AIC
mod_surv1 6 554.2907
mod_surv_null 1 564.5393show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(mod_surv1)
DHARMa::plotQQunif(mod_surv1)
Marginal \(R^2\) / Conditional \(R^2\)
show code
MuMIn::r.squaredGLMM(mod_surv1)Warning: the null model is only correct if all the variables it uses are identical
to those used in fitting the original model. R2m R2c
theoretical 0.1291593 0.4869551
delta 0.1160836 0.4376569show code
0.129/ 0.487[1] 0.2648871Which it is 0.265
show code
sjPlot::tab_model(mod_surv1)| cbind(n,24-n) | |||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.10 | 0.79 – 1.52 | 0.580 |
| tratamento [F25] | 0.80 | 0.62 – 1.03 | 0.090 |
| tratamento [F50] | 0.91 | 0.70 – 1.18 | 0.472 |
| tratamento [N25] | 0.58 | 0.45 – 0.75 | <0.001 |
| tratamento [N50] | 0.84 | 0.65 – 1.09 | 0.192 |
| Random Effects | |||
| σ2 | 3.29 | ||
| τ00 bloco | 0.10 | ||
| ICC | 0.03 | ||
| N bloco | 5 | ||
| Observations | 100 | ||
| Marginal R2 / Conditional R2 | 0.010 / 0.038 | ||
show code
broom.mixed::tidy(mod_surv1) |> dplyr::mutate(dplyr::across(4:7, round, 3)) Registered S3 method overwritten by 'broom':
method from
nobs.multinom MuMInWarning: There was 1 warning in `dplyr::mutate()`.
ℹ In argument: `dplyr::across(4:7, round, 3)`.
Caused by warning:
! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
Supply arguments directly to `.fns` through an anonymous function instead.
# Previously
across(a:b, mean, na.rm = TRUE)
# Now
across(a:b, \(x) mean(x, na.rm = TRUE))# A tibble: 6 × 7
effect group term estimate std.error statistic p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.092 0.166 0.553 0.58
2 fixed <NA> tratamentoF25 -0.222 0.131 -1.70 0.09
3 fixed <NA> tratamentoF50 -0.094 0.131 -0.719 0.472
4 fixed <NA> tratamentoN25 -0.544 0.132 -4.12 0
5 fixed <NA> tratamentoN50 -0.171 0.131 -1.31 0.192
6 ran_pars bloco sd__(Intercept) 0.309 NA NA NA show code
summary(mod_surv1)Generalized linear mixed model fit by maximum likelihood (Laplace
Approximation) [glmerMod]
Family: binomial ( logit )
Formula: cbind(n, 24 - n) ~ tratamento + (1 | bloco)
Data: survival_planted2024
AIC BIC logLik -2*log(L) df.resid
554.3 569.9 -271.1 542.3 94
Scaled residuals:
Min 1Q Median 3Q Max
-3.0363 -0.8992 -0.1404 1.0640 3.6310
Random effects:
Groups Name Variance Std.Dev.
bloco (Intercept) 0.0956 0.3092
Number of obs: 100, groups: bloco, 5
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.09194 0.16634 0.553 0.5805
tratamentoF25 -0.22197 0.13080 -1.697 0.0897 .
tratamentoF50 -0.09389 0.13068 -0.719 0.4724
tratamentoN25 -0.54411 0.13224 -4.115 3.88e-05 ***
tratamentoN50 -0.17069 0.13072 -1.306 0.1916
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) trtF25 trtF50 trtN25
tratamntF25 -0.393
tratamntF50 -0.393 0.500
tratamntN25 -0.389 0.494 0.495
tratamntN50 -0.393 0.500 0.500 0.495show code
exp(-0.544)[1] 0.5804219show code
(52.3/(100-52.3))/ (39.2/(100-39.2))[1] 1.700595Planted trees in treatment F25 have 19.91% less likely to survive (OR=exp(-0.222)) than in treatment 100% (\(p<0.1\))
Planted trees in treatment N25 have 41.96% less likely to survive (OR=exp(-0.544)) than in treatment 100% (\(p<0.001\))
show code
round((1 - exp(-0.544))*100,2)[1] 41.96show code
emmeans::emmeans(mod_surv1, ~ tratamento, type = "response") |>
multcomp::cld(details = TRUE, adjust="sidak")$emmeans
tratamento prob SE df asymp.LCL asymp.UCL .group
N25 0.389 0.0398 Inf 0.293 0.495 1
F25 0.468 0.0414 Inf 0.364 0.574 12
N50 0.480 0.0415 Inf 0.376 0.586 2
F50 0.500 0.0416 Inf 0.394 0.605 2
100 0.523 0.0415 Inf 0.417 0.627 2
Confidence level used: 0.95
Conf-level adjustment: sidak method for 5 estimates
Intervals are back-transformed from the logit scale
P value adjustment: sidak method for 10 tests
Tests are performed on the log odds ratio scale
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same.
$comparisons
contrast odds.ratio SE df null z.ratio p.value
F25 / N25 1.38 0.183 Inf 1 2.436 0.1390
N50 / N25 1.45 0.192 Inf 1 2.825 0.0463
N50 / F25 1.05 0.138 Inf 1 0.392 1.0000
F50 / N25 1.57 0.207 Inf 1 3.407 0.0066
F50 / F25 1.14 0.149 Inf 1 0.980 0.9810
F50 / N50 1.08 0.141 Inf 1 0.588 0.9997
100 / N25 1.72 0.228 Inf 1 4.115 0.0004
100 / F25 1.25 0.163 Inf 1 1.697 0.6093
100 / N50 1.19 0.155 Inf 1 1.306 0.8809
100 / F50 1.10 0.144 Inf 1 0.719 0.9983
P value adjustment: sidak method for 10 tests
Tests are performed on the log odds ratio scale (iii) Table S2 - Survival
show code
survival_planted2024 |>
dplyr::group_by( tratamento)|>
dplyr::summarise(Mean = round(mean((n/24)*100),1),
Max = round(max((n/24)*100),1),
Min = round(min((n/24)*100),1) ) |>
dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
dplyr::select(-c(Mean, Max, Min)) # A tibble: 5 × 2
tratamento values
<chr> <chr>
1 100 52.3 (29.2 - 79.2)
2 F25 46.9 (12.5 - 70.8)
3 F50 50 (29.2 - 75)
4 N25 39.2 (8.3 - 79.2)
5 N50 48.1 (20.8 - 75) B) Grass cover
(i) Organizing data
How planting scheme affected grass cover inside and outside planted areas?
Excluding treatment 100:
show code
ordem <- c("N25", "F25","N50", "F50", "100")
grass_cover <- readxl::read_xlsx(
paste0(path, "/Itatinga_Nuc_grasscover.xlsx"))|>
dplyr::mutate(
grass_perc = dplyr::if_else(gramineas_perc>=1, 0.9999, gramineas_perc),
tratamento = ordered(tratamento,
levels = ordem,
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%"))) |>
droplevels()
grass_cover2024 <- grass_cover |> dplyr::filter(tratamento!="Planting 100%") |> dplyr::filter(ano==2024)
grass_cover2024 |>
with(table(bloco, parcela, by=tratamento)), , by = Nuclei 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Nuclei 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Planting 100%
parcela
bloco 1 2 3 4
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0Too many zeros? Zero-Inflated model? There is few zeros. Only binomial fit
show code
grass_cover2024 |>
ggplot(aes(x=grass_perc))+
geom_histogram()+
facet_grid(plantio~tratamento)+theme_classic()`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
(ii) Analysis
show code
grass_cover2024 |>
with(table(bloco, parcela, by=tratamento)), , by = Nuclei 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 25%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Nuclei 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Strips 50%
parcela
bloco 1 2 3 4
1 2 2 2 2
2 2 2 2 2
3 2 2 2 2
4 2 2 2 2
5 2 2 2 2
, , by = Planting 100%
parcela
bloco 1 2 3 4
1 0 0 0 0
2 0 0 0 0
3 0 0 0 0
4 0 0 0 0
5 0 0 0 0show code
grass_cover2024 <-grass_cover2024 |>
dplyr::mutate(Tratamento=factor(tratamento, ordered = FALSE))
modelo_gc24a <- glmmTMB::glmmTMB(grass_perc ~ Tratamento*plantio + (1 | bloco),
family = glmmTMB::beta_family(), data = grass_cover2024)
modelo_gc24b <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ tratamento+plantio + (1 | bloco),
ziformula = ~ 1, family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + plantio + (1 | bloco), :
use of the 'data' argument is recommendedshow code
modelo_gc24c <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ plantio + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommendedshow code
modelo_gc24d <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ tratamento + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommendedshow code
modelo_gc24null <- grass_cover2024 |>
dplyr::mutate(tratamento=factor(tratamento,ordered = FALSE))|>
with(glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(grass_perc ~ 0 + (1 | bloco), family =
glmmTMB::beta_family()): use of the 'data' argument is recommendedshow code
AIC(modelo_gc24a, modelo_gc24b, modelo_gc24c, modelo_gc24d, modelo_gc24null) df AIC
modelo_gc24a 10 -61.25388
modelo_gc24b 8 -53.84400
modelo_gc24c 4 -51.58877
modelo_gc24d 6 -21.66566
modelo_gc24null 2 -10.55579Analysis of residuals
show code
par(mfrow=c(1,2))
DHARMa::plotQQunif(modelo_gc24a)
DHARMa::plotResiduals(modelo_gc24a)
Is there any difference between treatments considering inside/outside?
show code
resultado_24 <- emmeans::emmeans(modelo_gc24a, ~ Tratamento|plantio) |>
multcomp::cld( details = TRUE)
resultado_24$emmeansplantio = CP:
Tratamento emmean SE df asymp.LCL asymp.UCL .group
Strips 50% -0.727 0.228 Inf -1.173 -0.280 1
Strips 25% -0.300 0.224 Inf -0.739 0.138 12
Nuclei 25% 0.138 0.223 Inf -0.299 0.575 23
Nuclei 50% 0.708 0.228 Inf 0.262 1.154 3
plantio = SP:
Tratamento emmean SE df asymp.LCL asymp.UCL .group
Strips 50% 0.937 0.231 Inf 0.485 1.389 1
Nuclei 50% 0.945 0.231 Inf 0.493 1.397 1
Nuclei 25% 1.036 0.232 Inf 0.581 1.490 1
Strips 25% 1.166 0.234 Inf 0.708 1.625 1
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Results are given on the log odds ratio (not the response) scale.
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: yes.
Is there any difference between inside/outside considering treatments?
show code
resultado_24a <- emmeans::emmeans(modelo_gc24a, ~ plantio|Tratamento) |>
multcomp::cld(details = TRUE)
resultado_24a$emmeansTratamento = Nuclei 25%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP 0.138 0.223 Inf -0.299 0.575 1
SP 1.036 0.232 Inf 0.581 1.490 2
Tratamento = Strips 25%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP -0.300 0.224 Inf -0.739 0.138 1
SP 1.166 0.234 Inf 0.708 1.625 2
Tratamento = Nuclei 50%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP 0.708 0.228 Inf 0.262 1.154 1
SP 0.945 0.231 Inf 0.493 1.397 1
Tratamento = Strips 50%:
plantio emmean SE df asymp.LCL asymp.UCL .group
CP -0.727 0.228 Inf -1.173 -0.280 1
SP 0.937 0.231 Inf 0.485 1.389 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: yes.
show code
broom.mixed::tidy(modelo_gc24a) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))|> knitr::kable()| effect | component | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 0.138 | 0.223 | 0.620 | 0.535 |
| fixed | cond | NA | TratamentoStrips 25% | -0.439 | 0.316 | -1.388 | 0.165 |
| fixed | cond | NA | TratamentoNuclei 50% | 0.569 | 0.318 | 1.790 | 0.073 |
| fixed | cond | NA | TratamentoStrips 50% | -0.865 | 0.319 | -2.710 | 0.007 |
| fixed | cond | NA | plantioSP | 0.897 | 0.321 | 2.794 | 0.005 |
| fixed | cond | NA | TratamentoStrips 25%:plantioSP | 0.569 | 0.453 | 1.256 | 0.209 |
| fixed | cond | NA | TratamentoNuclei 50%:plantioSP | -0.660 | 0.454 | -1.454 | 0.146 |
| fixed | cond | NA | TratamentoStrips 50%:plantioSP | 0.766 | 0.454 | 1.688 | 0.091 |
| ran_pars | cond | bloco | sd__(Intercept) | 0.000 | NA | NA | NA |
show code
sjPlot::tab_model(modelo_gc24a)| grass perc | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 1.15 | 0.74 – 1.78 | 0.535 |
| Tratamento [Strips 25%] | 0.64 | 0.35 – 1.20 | 0.165 |
| Tratamento [Nuclei 50%] | 1.77 | 0.95 – 3.30 | 0.073 |
| Tratamento [Strips 50%] | 0.42 | 0.23 – 0.79 | 0.007 |
| plantio [SP] | 2.45 | 1.31 – 4.60 | 0.005 |
| Tratamento [Strips 25%] × plantio [SP] |
1.77 | 0.73 – 4.30 | 0.209 |
| Tratamento [Nuclei 50%] × plantio [SP] |
0.52 | 0.21 – 1.26 | 0.146 |
| Tratamento [Strips 50%] × plantio [SP] |
2.15 | 0.88 – 5.24 | 0.091 |
| Random Effects | |||
| σ2 | 0.16 | ||
| τ00 bloco | 0.00 | ||
| N bloco | 5 | ||
| Observations | 160 | ||
| Marginal R2 / Conditional R2 | 0.736 / NA | ||
show code
lme4::VarCorr(modelo_gc24a) # ou summary(modelo_gc24a)$varcor
Conditional model:
Groups Name Std.Dev.
bloco (Intercept) 1.0464e-05show code
MuMIn::r.squaredGLMM(modelo_gc24a)Warning: Can't compute random effect variances. Some variance components equal
zero. Your model may suffer from singularity (see `?lme4::isSingular`
and `?performance::check_singularity`).
Decrease the `tolerance` level to force the calculation of random effect
variances, or impose priors on your random effects parameters (using
packages like `brms` or `glmmTMB`). R2m R2c
[1,] 0.7355586 0.7355586show code
m_fix <- update(modelo_gc24a, . ~ . - (1|bloco))
anova(modelo_gc24a, m_fix) Data: grass_cover2024
Models:
m_fix: grass_perc ~ Tratamento + plantio + Tratamento:plantio, zi=~0, disp=~1
modelo_gc24a: grass_perc ~ Tratamento * plantio + (1 | bloco), zi=~0, disp=~1
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
m_fix 9 -63.254 -35.577 40.627 -81.254
modelo_gc24a 10 -61.254 -30.502 40.627 -81.254 0 1 1show code
performance::check_singularity(modelo_gc24a)[1] TRUE(iii) Figure - Grass cover
Small letters: Comparison between treatments, only inside planted areas.
Capital letters: Comparison between treatments, only outside planted areas.
show code
letras <- data.frame(
tratamento = rep(levels(as.factor(grass_cover$tratamento)),each=2),
plantio = levels(as.factor(grass_cover$plantio)),
label = c("BC ", " a",
"AB ", " a",
"C ", " a",
"A ", " a",
"A",""),
y = c(0.6, 0.82,
0.5, 0.8,
0.65, 0.8,
0.35, 0.8,
0.35, .1)
) |>
dplyr::mutate(
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
fig_grass <- grass_cover |>
dplyr::filter(ano==2024) |>
dplyr::group_by(tratamento, plantio) |>
dplyr::summarise(grass_cover_mean = mean(gramineas_perc),
sd = sd(gramineas_perc),
n = dplyr::n(),
se = sd/sqrt(n),
.groups="drop") |>
dplyr::mutate(
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
) |>
ggplot(aes(x = tratamento, y = grass_cover_mean, fill = plantio)) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = grass_cover_mean - se, ymax = grass_cover_mean + se),
width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('Planted', 'Unplanted'),
values = c("#92D050", "white")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size=4, vjust = -0.5) +
labs(fill = "",
y = "Grass cover (%)",
x = "") +
scale_y_continuous(expand = c(0,0),
labels = scales::percent_format(accuracy = 1),
limits = c(0,1)) +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=11),
axis.text.y = element_text(color="black", size=11),
axis.title.y = element_text(color="black", size=13),
legend.position = c(0.84, 0.99),
legend.text = element_text(size = 11, margin = margin(0,0,0, l=1, "mm")),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank())Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.show code
fig_grass
show code
# ggsave("fig_grass.png", plot=fig_grass, width=90, height = 90 ,unit="mm")(iv) Table S2 - Grass cover
show code
grass_cover |>
dplyr::group_by(ano, tratamento, plantio)|>
dplyr::summarise(Mean = mean(gramineas_perc)*100,
Max = max(gramineas_perc)*100,
Min = min(gramineas_perc)*100) |>
dplyr::mutate(values = paste0(Mean," (",Min," - ",Max,")")) |>
dplyr::select(-c(Mean, Max, Min)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)`summarise()` has grouped output by 'ano', 'tratamento'. You can override using
the `.groups` argument.# A tibble: 4 × 7
# Groups: ano [2]
ano plantio `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 CP 29.2 (0 - 85) 16.9 (0 - 100) 34.9 (0 - 100) 12.35 (0 - …
2 2022 SP 65.25 (20 - 100) 65.175 (23 - 100) 51.6 (0 - 100) 60.35 (0 - …
3 2024 CP 54.07 (8.2 - 95) 39.91 (4 - 95) 56.3485 (2.64… 29.9795 (5 …
4 2024 SP 76.835 (39 - 97.3) 75.45 (29 - 100) 73.627 (26 - … 77.9725 (49…
# ℹ 1 more variable: `Planting 100%` <chr>C) Recruitment
Only recruits (smaller and larger), not including planted individuals
show code
species_table <- readxl::read_xlsx("itatinga_Nuc_regen_species_abundance.xlsx")
species_table|>
dplyr::select(Treatment, Block, Plot, Planted)|>
unique()|>
with(table(Block, Plot, by=Treatment)), , by = Full_Planted
Plot
Block P1 P2 P3 P4
B1 1 1 1 1
B2 1 1 1 1
B3 1 1 1 1
B4 1 1 1 1
B5 1 1 1 1
, , by = Nuclei_25
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Nuclei_50
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Strip_25
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2
, , by = Strip_50
Plot
Block P1 P2 P3 P4
B1 2 2 2 2
B2 2 2 2 2
B3 2 2 2 2
B4 2 2 2 2
B5 2 2 2 2Viewing top species
show code
total_smaller <- sum(species_table$n_sapling)
total_larger <- sum(species_table$n_tree)
species_table |>
dplyr::group_by(scientificName)|>
dplyr::summarise(
sum_large = sum(n_tree, na.rm=T)/total_larger,
sum_small = sum(n_sapling, na.rm=T)/total_smaller,
)|>
dplyr::arrange(desc(sum_small)) |>
head()# A tibble: 6 × 3
scientificName sum_large sum_small
<chr> <dbl> <dbl>
1 Baccharis dracunculifolia 0.0740 0.489
2 Vernonanthura polyanthes 0 0.313
3 Solanum sisymbriifolium 0.00157 0.0388
4 Croton floribundus 0.00157 0.0384
5 Senegalia polyphylla 0.00472 0.0305
6 Solanum granuloso-leprosum 0.0205 0.0233(i) Considering Baccharis + Vernonanthura
Only small regenerants
show code
dados_list <- species_table |>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 2000, 100),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")Figure species richness
Figure Recruits’ species richness
For Hill q = 0
Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)
show code
dados_list <- species_table |>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 2000, 100),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")
curvas_df0 <- estimativas |>
dplyr::filter(Order.q==0)|>
dplyr::filter(!is.na(SC.LCL))|>
dplyr::mutate(
Percent = ifelse(Percent=="Planted", "100", Percent),
Percent = ordered(Percent, levels = c("25", "50", "100")),
Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
)
fig_recruit_rich <- curvas_df0 |>
ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
geom_point(
data = dplyr::filter(curvas_df0, Method == "Observed"),
size = 2.5
) +
facet_grid(Percent~.,
labeller = as_labeller(
c("25" = "(a) Planting 25%",
"50" = "(b) Planting 50%",
"100" = "(c) Planting 100%")
)) +
scale_color_manual(values=c("red", "blue", "darkgreen"))+
scale_fill_manual(values=c("red", "blue", "darkgreen"))+
labs(
x = "Number of sampling units (regenerants)",
y = "Species diversity (Hill Index, q=0)"
) +
ggthemes::theme_clean()+
scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 3500, 500)) +
scale_y_continuous(expand = c(0, 0))+
theme(
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color = "black", hjust = 1, size = 10),
axis.text.y = element_text(color = "black", size = 9),
axis.title.y = element_text(color = "black", size = 11),
strip.text = element_text(color = "black", hjust = 0.5, size = 10),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(5, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0, 0, 0, 0),
legend.position = c(0.25, 0.25),
panel.border = element_rect(colour = "black", fill = NA),
panel.spacing = unit(0, "lines")
)
fig_recruit_rich
show code
# ggsave("fig_richness_all_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")Table S2 - Recruitment (all)
Estimating Density of recruits (\(ind.ha^{-1}\)) and considering ruderal species (Baccharis and Vernonanthura)
show code
density_recruits_n_m2 <- species_table |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ha <- density_recruits_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_P_UP <- dens_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_regenerant_all <- dens_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
- Only in planted areas
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 3048.6 (2274.0 - 4705.0) 4871.5 (3490.0 - 7101.0) 6107.6 (2… 4816.0 … 6024.3 …- Only in unplanted areas
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 3961.8 (2344.0 - 8247.0) 3461.8 (2031.0 - … 3392.4 … 3871.5 …- In total (planted + unplanted areas)
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 3048.6 (2274.0 - 4705.0) 2094.6 (872.0 - 6185.0) 2392.4 (10… 1874.1 … 2474.0 …Larger recruits
All treatments and inside/outside areas showed similar values
- Only in planted areas
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 222.2 (104.0 - 486.0) 229.2 (0.0 - 608.0) 163.2 (35.0 - 295… 152.8 (… 225.7 (…- Only in unplanted areas
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 364.6 (17.0 - 851.0) 163.2 (122.0 - 208.0) 302.1 (1… 381.9 (…- In total (planted + unplanted areas)
show code
result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 222.2 (104.0 - 486.0) 165.4 (0.0 - 638.0) 81.6 (17.0 - 148.… 132.4 (… 151.9 (…Analysis for smaller recruits abundance
Data:
show code
small_df <- species_table |>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_sapling), .groups="drop")
small_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()
Models:
show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
# family = poisson(), # Poisson
family = glmmTMB::nbinom2(), # NB
data = small_df)
# NBinom 1782.688 (null) | 1734.985 (Treatment*Planted) | 1756.665 (only Treament) | 1743.097 (only Planted)
# Poisson 4032.177 | 3373.854 | 3813.769 | 3759.763
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 2 1782.688
modelo_reg1 11 1734.985
modelo_reg2 7 1756.665
modelo_reg3 4 1743.097non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:
show code
performance::check_singularity(modelo_reg1)[1] FALSEshow code
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)
show code
car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 23.4883 4 0.0001011 ***
Planted 30.8344 1 2.81e-08 ***
Treatment:Planted 1.8521 3 0.6036618
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Is there any difference between treatments considering only inside or outside?
show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust="bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.76 0.175 Inf 3.31 4.21 1
Strip_25 4.24 0.173 Inf 3.79 4.68 2
Nuclei_25 4.24 0.174 Inf 3.80 4.69 2
Nuclei_50 4.40 0.173 Inf 3.96 4.85 2
Strip_50 4.42 0.173 Inf 3.98 4.87 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 3.82 0.175 Inf 3.38 4.25 1
Nuclei_50 3.82 0.175 Inf 3.39 4.26 1
Nuclei_25 3.96 0.174 Inf 3.52 4.39 1
Strip_50 3.98 0.174 Inf 3.55 4.42 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: Yes (inside) and no (outside).
Is there any difference between treatments not considering inside/outside?
show code
emmeans::emmeans(modelo_reg2, ~ Treatment) |>
multcomp::cld(adjust="bonferroni") Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.76 0.177 Inf 3.30 4.22 1
Strip_25 4.05 0.155 Inf 3.65 4.45 12
Nuclei_25 4.11 0.155 Inf 3.71 4.51 12
Nuclei_50 4.16 0.155 Inf 3.76 4.56 12
Strip_50 4.23 0.154 Inf 3.83 4.63 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for 5 estimates
P value adjustment: bonferroni method for 10 tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Is there any difference between inside/outside considering treatments?
show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust="bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 3.76 0.175 Inf 3.42 4.10 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.96 0.174 Inf 3.57 4.35 1
Inside 4.24 0.174 Inf 3.85 4.63 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.82 0.175 Inf 3.43 4.21 1
Inside 4.40 0.173 Inf 4.02 4.79 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.82 0.175 Inf 3.42 4.21 1
Inside 4.24 0.173 Inf 3.85 4.63 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.98 0.174 Inf 3.59 4.37 1
Inside 4.42 0.173 Inf 4.03 4.81 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: Yes.
Table S3
show code
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.| n | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 42.99 | 30.52 – 60.55 | <0.001 |
| Treatment [Nuclei_25] | 1.62 | 1.19 – 2.20 | 0.002 |
| Treatment [Nuclei_50] | 1.90 | 1.40 – 2.59 | <0.001 |
| Treatment [Strip_25] | 1.61 | 1.19 – 2.19 | 0.002 |
| Treatment [Strip_50] | 1.94 | 1.43 – 2.63 | <0.001 |
| Planted [Outside] | 0.64 | 0.48 – 0.87 | 0.005 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.17 | 0.76 – 1.80 | 0.478 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.87 | 0.56 – 1.33 | 0.515 |
| Treatment [Strip_25] × Planted [Outside] |
1.02 | 0.66 – 1.57 | 0.937 |
| Random Effects | |||
| σ2 | 0.21 | ||
| τ00 Block | 0.09 | ||
| ICC | 0.30 | ||
| N Block | 5 | ||
| Observations | 180 | ||
| Marginal R2 / Conditional R2 | 0.164 / 0.412 | ||
show code
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))# A tibble: 11 × 6
effect term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed (Intercept) 3.76 0.175 21.5 0
2 fixed TreatmentNuclei_25 0.481 0.156 3.08 0.002
3 fixed TreatmentNuclei_50 0.644 0.156 4.13 0
4 fixed TreatmentStrip_25 0.476 0.156 3.05 0.002
5 fixed TreatmentStrip_50 0.662 0.156 4.25 0
6 fixed PlantedOutside -0.439 0.155 -2.83 0.005
7 fixed TreatmentNuclei_25:PlantedOuts… 0.156 0.219 0.709 0.478
8 fixed TreatmentNuclei_50:PlantedOuts… -0.143 0.22 -0.652 0.515
9 fixed TreatmentStrip_25:PlantedOutsi… 0.017 0.22 0.079 0.937
10 fixed TreatmentStrip_50:PlantedOutsi… NA NA NA NA
11 ran_pars sd__(Intercept) 0.301 NA NA NA Analysis for larger recruits abundance
Data:
show code
large_regenerants_2024 <- species_table |>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_tree), .groups="drop")Plotting
show code
large_regenerants_2024 |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()
Models
show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 885.72 | 869.63 | 868.93 | 864.21
# ZIP 1058.14 | 1015.06 | 1025.38 | 1024.49
# Hurdle 880.75 | 872.21 | 868.74 | 975.88
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 3 885.7221
modelo_reg1 12 869.6264
modelo_reg2 8 868.9342
modelo_reg3 5 864.2072show code
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)
show code
MuMIn::r.squaredGLMM(modelo_reg3)Warning: the null model is only correct if all the variables it uses are identical
to those used in fitting the original model. R2m R2c
delta 0.03380077 0.03380078
lognormal 0.05135140 0.05135140
trigamma 0.01854201 0.01854201show code
# 0.03253385/0.036655970.0328/ 0.0338 variance explained by fixed effects and by the entire model (around 0.9999997 % is explained by fixed effects)
show code
performance::r2(modelo_reg1)Warning: Can't compute random effect variances. Some variance components equal
zero. Your model may suffer from singularity (see `?lme4::isSingular`
and `?performance::check_singularity`).
Decrease the `tolerance` level to force the calculation of random effect
variances, or impose priors on your random effects parameters (using
packages like `brms` or `glmmTMB`).[1] NAshow code
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
dplyr::select(-c(component, group))|>
knitr::kable()| effect | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| fixed | (Intercept) | 1.175 | 0.251 | 4.690 | 0.000 |
| fixed | TreatmentNuclei_25 | 0.064 | 0.369 | 0.173 | 0.863 |
| fixed | TreatmentNuclei_50 | -0.286 | 0.369 | -0.776 | 0.438 |
| fixed | TreatmentStrip_25 | -0.353 | 0.369 | -0.957 | 0.339 |
| fixed | TreatmentStrip_50 | 0.028 | 0.354 | 0.080 | 0.936 |
| fixed | PlantedOutside | 0.510 | 0.346 | 1.476 | 0.140 |
| fixed | TreatmentNuclei_25:PlantedOutside | -0.045 | 0.492 | -0.092 | 0.927 |
| fixed | TreatmentNuclei_50:PlantedOutside | -0.515 | 0.507 | -1.018 | 0.309 |
| fixed | TreatmentStrip_25:PlantedOutside | 0.169 | 0.499 | 0.339 | 0.734 |
| fixed | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | (Intercept) | -3.510 | 3.051 | -1.150 | 0.250 |
| ran_pars | sd__(Intercept) | 0.000 | NA | NA | NA |
show code
car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 7.3216 4 0.11984
Planted 5.5546 1 0.01843 *
Treatment:Planted 1.9111 3 0.59106
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Is there any difference between treatments considering inside/outside?
show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 0.822 0.284 Inf 0.0899 1.55 1
Nuclei_50 0.889 0.284 Inf 0.1589 1.62 1
Full_Planted 1.175 0.251 Inf 0.5298 1.82 1
Strip_50 1.203 0.260 Inf 0.5327 1.87 1
Nuclei_25 1.239 0.287 Inf 0.4988 1.98 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 0.884 0.277 Inf 0.1920 1.58 1
Strip_25 1.501 0.259 Inf 0.8532 2.15 1
Nuclei_25 1.704 0.276 Inf 1.0145 2.39 1
Strip_50 1.714 0.235 Inf 1.1262 2.30 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Is there any difference between inside/outside considering treatments?
show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.175 0.251 Inf 0.684 1.67 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.239 0.287 Inf 0.595 1.88 1
Outside 1.704 0.276 Inf 1.085 2.32 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 0.884 0.277 Inf 0.263 1.51 1
Inside 0.889 0.284 Inf 0.254 1.52 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.822 0.284 Inf 0.185 1.46 1
Outside 1.501 0.259 Inf 0.920 2.08 1
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.203 0.260 Inf 0.620 1.79 1
Outside 1.714 0.235 Inf 1.186 2.24 1
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Propotion of smaller recruits from planted species
- With ruderal species
Getting proportions of planted/unplanted species
show code
prop_regenerant_full <- species_table |>
dplyr::mutate(
planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
dplyr::filter(n_sapling>0)|>
dplyr::group_by(Treatment, Planted, Block, Plot) |>
dplyr::summarise(
n_species = dplyr::n(),
n_ind = sum(n_sapling),
n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
n_planted_spp = sum(planted_tree == "Yes"),
n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
.groups = "drop") |>
dplyr::mutate(
prop_planted_ind = n_planted_ind / n_ind,
prop_planted_spp = n_planted_spp / n_species
)
prop_regenerant_full |>
dplyr::group_by(Treatment, Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_min = min(prop_planted_ind, na.rm=T),
prop_planted_max = max(prop_planted_ind, na.rm=T),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
.groups = "drop"
) |>
dplyr::mutate(
prop_planted = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_mean,3)*100,
round(prop_planted_min,3)*100,
round(prop_planted_max,3)*100),
prop_planted_spp = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_spp_mean,3)*100,
round(prop_planted_spp_min,3)*100,
round(prop_planted_spp_max,3)*100)
) |>
dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = c(prop_planted, prop_planted_spp)
) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")# A tibble: 2 × 11
Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
<chr> <chr> <chr> <chr>
1 Inside 10.5 (0.0 - 43.3) 3.9 (0.0 - 14.8) 7.7 (0.0 - 44.4)
2 Outside <NA> 10.4 (0.0 - 91.2) 10.0 (0.0 - 69.2)
# ℹ abbreviated name: ¹prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
# prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
# prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
# prop_planted_spp_Strip_50 <chr>Statistical analysis
show code
mod_binomial_full <-lme4::glmer(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_full)fixed-effect model matrix is rank deficient so dropping 1 column / coefficientshow code
performance::check_overdispersion(mod_binomial_full) # no# Overdispersion test
dispersion ratio = 1.748
p-value = 0.272No overdispersion detected.show code
performance::check_model(mod_binomial_full)
show code
DHARMa::plotResiduals(mod_binomial_full)
show code
emmeans::emmeans(mod_binomial_full, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -3.35 0.292 Inf -4.10 -2.598 1
Strip_50 -2.63 0.271 Inf -3.32 -1.933 2
Nuclei_50 -2.59 0.269 Inf -3.28 -1.901 2
Strip_25 -2.40 0.269 Inf -3.09 -1.708 2
Full_Planted -2.30 0.279 Inf -3.02 -1.583 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 -2.79 0.285 Inf -3.50 -2.084 1
Nuclei_25 -2.48 0.278 Inf -3.17 -1.789 12
Nuclei_50 -2.30 0.276 Inf -2.99 -1.610 2
Strip_50 -1.56 0.265 Inf -2.22 -0.903 3
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
emmeans::emmeans(mod_binomial_full, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.30 0.279 Inf -2.85 -1.753 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -3.35 0.292 Inf -4.00 -2.695 1
Outside -2.48 0.278 Inf -3.10 -1.860 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.59 0.269 Inf -3.19 -1.991 1
Outside -2.30 0.276 Inf -2.92 -1.681 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -2.79 0.285 Inf -3.43 -2.156 1
Inside -2.40 0.269 Inf -3.00 -1.797 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside -2.63 0.271 Inf -3.23 -2.023 1
Outside -1.56 0.265 Inf -2.16 -0.971 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
broom.mixed::tidy(mod_binomial_full) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 10 × 7
effect group term estimate std.error statistic p.value
<chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -2.30 0.279 -8.24 0
2 fixed <NA> TreatmentNuclei_25 -1.05 0.188 -5.56 0
3 fixed <NA> TreatmentNuclei_50 -0.293 0.151 -1.95 0.051
4 fixed <NA> TreatmentStrip_25 -0.101 0.151 -0.665 0.506
5 fixed <NA> TreatmentStrip_50 -0.329 0.154 -2.14 0.033
6 fixed <NA> PlantedOutside 1.06 0.128 8.33 0
7 fixed <NA> TreatmentNuclei_25:Plant… -0.199 0.226 -0.88 0.379
8 fixed <NA> TreatmentNuclei_50:Plant… -0.77 0.193 -3.99 0
9 fixed <NA> TreatmentStrip_25:Plante… -1.46 0.205 -7.12 0
10 ran_pars Block sd__(Intercept) 0.564 NA NA NA show code
sjPlot::tab_model(mod_binomial_full)| cbind(n planted ind,n ind-n planted ind) |
|||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 0.10 | 0.06 – 0.17 | <0.001 |
| Treatment [Nuclei_25] | 0.35 | 0.24 – 0.51 | <0.001 |
| Treatment [Nuclei_50] | 0.75 | 0.56 – 1.00 | 0.051 |
| Treatment [Strip_25] | 0.90 | 0.67 – 1.22 | 0.506 |
| Treatment [Strip_50] | 0.72 | 0.53 – 0.97 | 0.033 |
| Planted [Outside] | 2.90 | 2.26 – 3.72 | <0.001 |
| Treatment [Nuclei_25] × Planted [Outside] |
0.82 | 0.53 – 1.28 | 0.379 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.46 | 0.32 – 0.68 | <0.001 |
| Treatment [Strip_25] × Planted [Outside] |
0.23 | 0.16 – 0.35 | <0.001 |
| Random Effects | |||
| σ2 | 0.05 | ||
| τ00 Block | 0.32 | ||
| ICC | 0.86 | ||
| N Block | 5 | ||
| Observations | 178 | ||
| Marginal R2 / Conditional R2 | 0.354 / 0.910 | ||
(ii) Not considering Baccharis + Vernonanthura
Figure species richness
Hill q = 0
Only small regenerants and considering ruderal species (Vernonanthura & Baccharis)
show code
dados_list <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Planted, scientificName) |>
dplyr::summarise(n_total = sum(n_sapling), .groups = "drop") |>
dplyr::group_by(Treatment, Planted) |>
dplyr::summarise(abundancias = list(n_total), .groups = "drop") |>
dplyr::mutate(nome_grupo = paste(Treatment, Planted, sep = "_")) %>%
{ setNames(.$abundancias, .$nome_grupo) }
res_inext <- iNEXT::iNEXT(
dados_list,
q = 0,
size=seq(0, 450, 10),
datatype = "abundance"
)
estimativas <- res_inext$iNextEst |>
dplyr::bind_rows() |>
tidyr::separate(Assemblage, into = c("Treatment", "Percent", "Planted"), sep = "_")
curvas_df0 <- estimativas |>
dplyr::filter(Order.q==0)|>
dplyr::filter(!is.na(SC.LCL))|>
dplyr::mutate(
Percent = ifelse(Percent=="Planted", "100", Percent),
Percent = ordered(Percent, levels = c("25", "50", "100")),
Planted = ifelse(Planted=="Inside", "Planted", "Unplanted"),
Treatment = ifelse(Treatment=="Full", "Planting 100%", Treatment),
Treatment = factor(Treatment, levels = c("Nuclei", "Strip", "Planting 100%"), ordered = TRUE)
)
fig_recruit_rich <- curvas_df0 |>
ggplot(aes(x = m, y = qD, color = Treatment, linetype = Planted)) +
geom_line(linewidth = 1) +
geom_ribbon(aes(ymin = qD.LCL, ymax = qD.UCL, fill = Treatment), alpha = 0.2, colour = NA) +
geom_point(
data = dplyr::filter(curvas_df0, Method == "Observed"),
size = 2.5
) +
facet_grid(Percent~.,
labeller = as_labeller(
c("25" = "(a) Planting 25%",
"50" = "(b) Planting 50%",
"100" = "(c) Planting 100%")
)) +
scale_color_manual(values=c("red", "blue", "darkgreen"))+
scale_fill_manual(values=c("red", "blue", "darkgreen"))+
labs(
x = "Number of sampling units (regenerants)",
y = "Species diversity (Hill Index, q=0)"
) +
ggthemes::theme_clean()+
scale_x_continuous(expand = c(0.001, 0), breaks = seq(0, 450, 50)) +
scale_y_continuous(expand = c(0, 0))+
theme(
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color = "black", hjust = 1, size = 10),
axis.text.y = element_text(color = "black", size = 9),
axis.title.y = element_text(color = "black", size = 11),
strip.text = element_text(color = "black", hjust = 0.5, size = 10),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(5, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0, 0, 0, 0),
legend.position = c(0.25, 0.25),
panel.border = element_rect(colour = "black", fill = NA),
panel.spacing = unit(0, "lines")
)Plotting and saving
show code
fig_recruit_rich
show code
ggsave("fig_richness_no_ruderal_small_recruits.png", plot=fig_recruit_rich, width=90, height = 210 ,unit="mm")Table S2 - Recruitment (no ruderal)
Estimating Density of recruits (\(ind.ha^{-1}\)) without Baccharis and Vernonanthura
show code
density_recruits_n_m2 <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ha <- density_recruits_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_P_UP <- dens_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_regenerant <- dens_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
- Only in planted areas (without ruderal species)
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 503.5 (139.0 - 1146.0) 482.6 (174.0 - 903.0) 961.8 (278.0 -… 722.2 (… 614.6 (…- Only in unplanted areas (without ruderal species)
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 1177.1 (469.0 - 1910.0) 822.9 (538.0 - 138… 1145.8 … 1381.9 …- In total (planted + unplanted areas)
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 503.5 (139.0 - 1146.0) 501.7 (43.0 - 1432.0) 446.2 (139.0 -… 577.7 (… 499.1 (…Larger recruits
- Only in planted areas (without ruderal species)
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
) |>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 208.3 (87.0 - 486.0) 211.8 (0.0 - 608.0) 159.7 (35.0 - 278.… 149.3 (… 222.2 (…- Only in unplanted areas
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 333.3 (0.0 - 851.0) 138.9 (104.0 - 208.0) 253.5 (35… 364.6 (…- In total (planted + unplanted areas)
show code
result_regenerant |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 208.3 (87.0 - 486.0) 151.5 (0.0 - 638.0) 74.7 (17.0 - 139.0) 126.4 (… 146.7 (…Analysis for smaller recruits abundance
Data:
show code
small_df <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_sapling), .groups="drop")
small_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()
Models:
show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
# family = poisson(),
family = glmmTMB::nbinom2(),
data = small_df)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
# family = poisson(), # Poisson
family = glmmTMB::nbinom2(), # NB
data = small_df)
# NBinom 1270.77 (null) | 1226.20 (Treatment*Planted) | 1250.26 (only Treament) | 1227.00 (only Planted)
# Poisson 4032.177 | 3373.854 | 3813.769 | 3759.763
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 2 1270.773
modelo_reg1 11 1226.199
modelo_reg2 7 1250.263
modelo_reg3 4 1227.001non-positive-definite Hessian matrix indicates that there are problems in estimating the model parameters. This may be related to redundancies in the random effects, perfect separation, few data for some combinations of factors, or problems in choosing the distribution. Therefore, let’s check the singularity:
show code
performance::check_singularity(modelo_reg1)[1] FALSEshow code
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)
show code
car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 4.7099 4 0.31838
Planted 24.8150 1 6.311e-07 ***
Treatment:Planted 10.7731 3 0.01302 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Is there any difference between treatments considering only inside or outside?
show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust="bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 1.86 0.241 Inf 1.24 2.48 1
Nuclei_25 1.87 0.241 Inf 1.25 2.50 1
Strip_50 2.09 0.238 Inf 1.48 2.70 1
Strip_25 2.27 0.236 Inf 1.66 2.88 1
Nuclei_50 2.42 0.234 Inf 1.82 3.02 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 2.43 0.234 Inf 1.84 3.01 1
Strip_25 2.73 0.232 Inf 2.15 3.31 1
Nuclei_25 2.79 0.231 Inf 2.21 3.36 1
Strip_50 2.99 0.230 Inf 2.41 3.56 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No
Is there any difference between inside/outside considering treatments?
show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust="bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.86 0.241 Inf 1.39 2.33 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.87 0.241 Inf 1.33 2.42 1
Outside 2.79 0.231 Inf 2.27 3.30 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.42 0.234 Inf 1.90 2.95 1
Outside 2.43 0.234 Inf 1.90 2.95 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.27 0.236 Inf 1.74 2.80 1
Outside 2.73 0.232 Inf 2.21 3.25 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 2.09 0.238 Inf 1.56 2.63 1
Outside 2.99 0.230 Inf 2.47 3.51 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: Yes.
show code
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.| n | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p |
| (Intercept) | 6.43 | 4.01 – 10.31 | <0.001 |
| Treatment [Nuclei_25] | 1.01 | 0.63 – 1.62 | 0.955 |
| Treatment [Nuclei_50] | 1.75 | 1.11 – 2.76 | 0.016 |
| Treatment [Strip_25] | 1.51 | 0.95 – 2.39 | 0.079 |
| Treatment [Strip_50] | 1.26 | 0.79 – 2.00 | 0.326 |
| Planted [Outside] | 2.45 | 1.57 – 3.82 | <0.001 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.01 | 0.54 – 1.91 | 0.964 |
| Treatment [Nuclei_50] × Planted [Outside] |
0.41 | 0.22 – 0.77 | 0.005 |
| Treatment [Strip_25] × Planted [Outside] |
0.65 | 0.35 – 1.21 | 0.171 |
| Random Effects | |||
| σ2 | 0.40 | ||
| τ00 Block | 0.15 | ||
| ICC | 0.27 | ||
| N Block | 5 | ||
| Observations | 180 | ||
| Marginal R2 / Conditional R2 | 0.207 / 0.419 | ||
show code
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |> dplyr::select(-c(group,component))# A tibble: 11 × 6
effect term estimate std.error statistic p.value
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed (Intercept) 1.86 0.241 7.72 0
2 fixed TreatmentNuclei_25 0.014 0.24 0.056 0.955
3 fixed TreatmentNuclei_50 0.56 0.232 2.41 0.016
4 fixed TreatmentStrip_25 0.411 0.234 1.75 0.079
5 fixed TreatmentStrip_50 0.232 0.236 0.982 0.326
6 fixed PlantedOutside 0.897 0.226 3.97 0
7 fixed TreatmentNuclei_25:PlantedOuts… 0.015 0.322 0.045 0.964
8 fixed TreatmentNuclei_50:PlantedOuts… -0.892 0.319 -2.80 0.005
9 fixed TreatmentStrip_25:PlantedOutsi… -0.437 0.319 -1.37 0.171
10 fixed TreatmentStrip_50:PlantedOutsi… NA NA NA NA
11 ran_pars sd__(Intercept) 0.383 NA NA NA Analysis for larger recruits abundance
Data:
show code
large_regenerants_2024 <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_tree), .groups="drop")Plotting
show code
large_regenerants_2024 |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()
Models
show code
modelo_reg1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
modelo_reg2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = large_regenerants_2024)
modelo_reg_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = large_regenerants_2024)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZINB 855.21 | 843.24 | 840.70 | 838.42
AIC(modelo_reg_null, modelo_reg1, modelo_reg2, modelo_reg3) df AIC
modelo_reg_null 3 855.2141
modelo_reg1 12 843.2460
modelo_reg2 8 840.7029
modelo_reg3 5 838.4200show code
performance::check_model(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotQQunif(modelo_reg1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotResiduals(modelo_reg3)
DHARMa::plotQQunif(modelo_reg3)
show code
car::Anova(modelo_reg1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 7.9554 4 0.09323 .
Planted 3.4787 1 0.06216 .
Treatment:Planted 2.1403 3 0.54381
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Is there any difference between treatments considering inside/outside?
show code
emmeans::emmeans(modelo_reg1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 0.920 0.292 Inf 0.170 1.67 1
Nuclei_50 0.971 0.282 Inf 0.245 1.70 1
Full_Planted 1.162 0.249 Inf 0.520 1.80 1
Strip_50 1.278 0.261 Inf 0.606 1.95 1
Nuclei_25 1.280 0.280 Inf 0.558 2.00 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 0.820 0.284 Inf 0.109 1.53 1
Strip_25 1.419 0.259 Inf 0.772 2.07 1
Strip_50 1.744 0.238 Inf 1.150 2.34 1
Nuclei_25 1.787 0.283 Inf 1.079 2.49 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Is there any difference between inside/outside considering treatments?
show code
emmeans::emmeans(modelo_reg1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.162 0.249 Inf 0.674 1.65 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.280 0.280 Inf 0.652 1.91 1
Outside 1.787 0.283 Inf 1.152 2.42 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 0.820 0.284 Inf 0.182 1.46 1
Inside 0.971 0.282 Inf 0.339 1.60 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.920 0.292 Inf 0.267 1.57 1
Outside 1.419 0.259 Inf 0.838 2.00 1
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 1.278 0.261 Inf 0.694 1.86 1
Outside 1.744 0.238 Inf 1.211 2.28 1
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Table S3
show code
sjPlot::tab_model(modelo_reg1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.| n | ||||
|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p | |
| Count Model | ||||
| (Intercept) | 3.20 | 1.96 – 5.21 | <0.001 | |
| Treatment [Nuclei_25] | 1.13 | 0.55 – 2.31 | 0.747 | |
| Treatment [Nuclei_50] | 0.83 | 0.40 – 1.70 | 0.605 | |
| Treatment [Strip_25] | 0.79 | 0.38 – 1.64 | 0.520 | |
| Treatment [Strip_50] | 1.12 | 0.56 – 2.25 | 0.743 | |
| Planted [Outside] | 1.59 | 0.81 – 3.14 | 0.178 | |
| Treatment [Nuclei_25] × Planted [Outside] |
1.04 | 0.38 – 2.82 | 0.937 | |
| Treatment [Nuclei_50] × Planted [Outside] |
0.54 | 0.20 – 1.49 | 0.233 | |
| Treatment [Strip_25] × Planted [Outside] |
1.03 | 0.38 – 2.80 | 0.949 | |
| (Intercept) | 3.53 | 2.03 – 9.52 | ||
| Zero-Inflated Model | ||||
| (Intercept) | 0.14 | 0.04 – 0.50 | 0.003 | |
| Random Effects | ||||
| σ2 | ||||
| τ00 Block | 0.00 | |||
| N Block | 5 | |||
| Observations | 180 | |||
show code
broom.mixed::tidy(modelo_reg1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
knitr::kable()| effect | component | group | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 1.162 | 0.249 | 4.663 | 0.000 |
| fixed | cond | NA | TreatmentNuclei_25 | 0.119 | 0.367 | 0.323 | 0.747 |
| fixed | cond | NA | TreatmentNuclei_50 | -0.191 | 0.369 | -0.517 | 0.605 |
| fixed | cond | NA | TreatmentStrip_25 | -0.241 | 0.376 | -0.643 | 0.520 |
| fixed | cond | NA | TreatmentStrip_50 | 0.116 | 0.355 | 0.328 | 0.743 |
| fixed | cond | NA | PlantedOutside | 0.466 | 0.346 | 1.347 | 0.178 |
| fixed | cond | NA | TreatmentNuclei_25:PlantedOutside | 0.040 | 0.508 | 0.080 | 0.937 |
| fixed | cond | NA | TreatmentNuclei_50:PlantedOutside | -0.617 | 0.518 | -1.192 | 0.233 |
| fixed | cond | NA | TreatmentStrip_25:PlantedOutside | 0.033 | 0.508 | 0.064 | 0.949 |
| fixed | cond | NA | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | zi | NA | (Intercept) | -1.962 | 0.652 | -3.007 | 0.003 |
| ran_pars | cond | Block | sd__(Intercept) | 0.000 | NA | NA | NA |
Propotion of smaller recruits from planted species
- Without ruderal species
Getting proportions of planted/unplanted species
show code
prop_regenerant_part <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
dplyr::mutate(
planted_tree = ifelse(`Planting group`=="NA", "No", "Yes"))|>
dplyr::filter(n_sapling>0)|>
dplyr::group_by(Treatment, Planted, Block, Plot) |>
dplyr::summarise(
n_species = dplyr::n(),
n_ind = sum(n_sapling),
n_planted_ind = sum(ifelse(planted_tree == "Yes", n_sapling, 0), na.rm = TRUE),
n_planted_ind = ifelse(is.na(n_planted_ind), 0, n_planted_ind),
n_planted_spp = sum(planted_tree == "Yes"),
n_planted_spp = ifelse(is.na(n_planted_spp), 0, n_planted_spp),
.groups = "drop") |>
dplyr::mutate(
prop_planted_ind = n_planted_ind / n_ind,
prop_planted_spp = n_planted_spp / n_species
)
prop_regenerant_part |>
dplyr::group_by(Treatment, Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_min = min(prop_planted_ind, na.rm=T),
prop_planted_max = max(prop_planted_ind, na.rm=T),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_min = min(prop_planted_spp, na.rm=T),
prop_planted_spp_max = max(prop_planted_spp, na.rm=T),
.groups = "drop"
) |>
dplyr::mutate(
prop_planted = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_mean,3)*100,
round(prop_planted_min,3)*100,
round(prop_planted_max,3)*100),
prop_planted_spp = sprintf("%.1f (%.1f - %.1f)",
round(prop_planted_spp_mean,3)*100,
round(prop_planted_spp_min,3)*100,
round(prop_planted_spp_max,3)*100)
) |>
dplyr::select(Treatment,Planted, prop_planted, prop_planted_spp) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = c(prop_planted, prop_planted_spp)
) # |> writexl::write_xlsx("resulta_prop_recruit.xlsx")# A tibble: 2 × 11
Planted prop_planted_Full_Plan…¹ prop_planted_Nuclei_25 prop_planted_Nuclei_50
<chr> <chr> <chr> <chr>
1 Inside 54.4 (0.0 - 100.0) 41.4 (0.0 - 100.0) 37.5 (0.0 - 100.0)
2 Outside <NA> 25.0 (0.0 - 91.2) 27.0 (0.0 - 92.3)
# ℹ abbreviated name: ¹prop_planted_Full_Planted
# ℹ 7 more variables: prop_planted_Strip_25 <chr>, prop_planted_Strip_50 <chr>,
# prop_planted_spp_Full_Planted <chr>, prop_planted_spp_Nuclei_25 <chr>,
# prop_planted_spp_Nuclei_50 <chr>, prop_planted_spp_Strip_25 <chr>,
# prop_planted_spp_Strip_50 <chr>show code
prop_regenerant_part |>
dplyr::group_by(Planted)|>
dplyr::summarise(
prop_planted_mean = mean(prop_planted_ind, na.rm=T),
prop_planted_se = (sd(prop_planted_ind, na.rm=T))/sqrt(dplyr::n()),
prop_planted_spp_mean = mean(prop_planted_spp, na.rm=T),
prop_planted_spp_se = (sd(prop_planted_spp, na.rm=T))/sqrt(dplyr::n()),
.groups = "drop"
)# A tibble: 2 × 5
Planted prop_planted_mean prop_planted_se prop_planted_spp_mean
<chr> <dbl> <dbl> <dbl>
1 Inside 0.524 0.0372 0.488
2 Outside 0.283 0.0309 0.243
# ℹ 1 more variable: prop_planted_spp_se <dbl>Statistical analysis for proportion of individuals:
show code
mod_binomial_part <-lme4::glmer(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_part)fixed-effect model matrix is rank deficient so dropping 1 column / coefficientshow code
performance::check_overdispersion(mod_binomial_part) # yes# Overdispersion test
dispersion ratio = 1.837
p-value = < 0.001Overdispersion detected.It should be Beta
show code
mod_beta_part <-glmmTMB::glmmTMB(
cbind(n_planted_ind, n_ind - n_planted_ind) ~ Treatment*Planted + (1 | Block),
family = glmmTMB::betabinomial(link = "logit"),
data = prop_regenerant_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
performance::check_model(mod_beta_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
show code
DHARMa::plotResiduals(mod_beta_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Pairwise:
show code
emmeans::emmeans(mod_beta_part, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differencesPlanted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_50 -0.429 0.311 Inf -1.228 0.370 1
Nuclei_25 -0.358 0.314 Inf -1.165 0.449 1
Full_Planted 0.243 0.326 Inf -0.595 1.082 1
Strip_25 0.447 0.298 Inf -0.319 1.212 1
Strip_50 0.490 0.292 Inf -0.260 1.240 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 -1.196 0.298 Inf -1.938 -0.455 1
Nuclei_25 -1.096 0.287 Inf -1.810 -0.382 1
Nuclei_50 -0.935 0.314 Inf -1.718 -0.153 1
Strip_50 -0.376 0.264 Inf -1.035 0.282 1
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
emmeans::emmeans(mod_beta_part, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak") # All differentTreatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.243 0.326 Inf -0.397 0.883 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.096 0.287 Inf -1.737 -0.455 1
Inside -0.358 0.314 Inf -1.061 0.345 1
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.935 0.314 Inf -1.638 -0.233 1
Inside -0.429 0.311 Inf -1.125 0.267 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.196 0.298 Inf -1.862 -0.530 1
Inside 0.447 0.298 Inf -0.220 1.113 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.376 0.264 Inf -0.968 0.215 1
Inside 0.490 0.292 Inf -0.163 1.143 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Statistical analysis for proportion of individuals:
show code
mod_binomial_prop <-lme4::glmer(
cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
family = binomial(link = "logit"),
data = prop_regenerant_part)fixed-effect model matrix is rank deficient so dropping 1 column / coefficientboundary (singular) fit: see help('isSingular')show code
performance::check_overdispersion(mod_binomial_prop) # underdispersion# Overdispersion test
dispersion ratio = 0.769
p-value = 0.032Underdispersion detected.It should be Beta
show code
mod_beta_prop <-glmmTMB::glmmTMB(
cbind(n_planted_spp, n_species - n_planted_spp) ~ Treatment*Planted + (1 | Block),
family = glmmTMB::betabinomial(link = "logit"),
data = prop_regenerant_part)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideWarning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; non-positive-definite Hessian matrix. See vignette('troubleshooting')Warning in finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old): Model convergence
problem; false convergence (8). See vignette('troubleshooting'),
help('diagnose')show code
performance::check_model(mod_beta_prop)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
show code
DHARMa::plotResiduals(mod_beta_prop)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside
Pairwise:
show code
emmeans::emmeans(mod_beta_prop, pairwise ~ Treatment | Planted) |> multcomp::cld(adjust="sidak") # no differencesPlanted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -0.582 0.286 Inf -1.318 0.154 1
Nuclei_50 -0.528 0.263 Inf -1.203 0.147 1
Strip_25 -0.210 0.246 Inf -0.841 0.421 1
Full_Planted 0.251 0.291 Inf -0.496 0.999 1
Strip_50 0.297 0.259 Inf -0.368 0.962 1
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Nuclei_25 -1.399 0.250 Inf -2.021 -0.777 1
Strip_25 -1.204 0.233 Inf -1.784 -0.624 1
Nuclei_50 -0.944 0.257 Inf -1.585 -0.304 1
Strip_50 -0.847 0.208 Inf -1.366 -0.329 1
Full_Planted nonEst NA NA NA NA
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
P value adjustment: sidak method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
emmeans::emmeans(mod_beta_prop, pairwise ~ Planted | Treatment) |> multcomp::cld(adjust="sidak") # All differentTreatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 0.251 0.291 Inf -0.319 0.8216 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.399 0.250 Inf -1.957 -0.8403 1
Inside -0.582 0.286 Inf -1.222 0.0587 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.944 0.257 Inf -1.520 -0.3693 1
Inside -0.528 0.263 Inf -1.116 0.0599 1
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -1.204 0.233 Inf -1.724 -0.6835 1
Inside -0.210 0.246 Inf -0.759 0.3398 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside -0.847 0.208 Inf -1.313 -0.3820 1
Inside 0.297 0.259 Inf -0.282 0.8763 2
Results are given on the logit (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: sidak method for varying numbers of estimates
Results are given on the log odds ratio (not the response) scale.
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Table S3
show code
broom.mixed::tidy(mod_beta_part) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 11 × 8
effect component group term estimate std.error statistic p.value
<chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 fixed cond <NA> (Intercept) 0.243 0.326 0.745 0.456
2 fixed cond <NA> TreatmentNucle… -0.601 0.454 -1.32 0.185
3 fixed cond <NA> TreatmentNucle… -0.672 0.451 -1.49 0.136
4 fixed cond <NA> TreatmentStrip… 0.203 0.442 0.461 0.645
5 fixed cond <NA> TreatmentStrip… 0.247 0.437 0.564 0.572
6 fixed cond <NA> PlantedOutside -0.867 0.394 -2.20 0.028
7 fixed cond <NA> TreatmentNucle… 0.129 0.577 0.223 0.824
8 fixed cond <NA> TreatmentNucle… 0.36 0.592 0.608 0.543
9 fixed cond <NA> TreatmentStrip… -0.776 0.576 -1.35 0.178
10 fixed cond <NA> TreatmentStrip… NA NA NA NA
11 ran_pars cond Block sd__(Intercept) 0 NA NA NA show code
sjPlot::tab_model(mod_beta_part)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.| cbind(n planted ind,n ind-n planted ind) |
|||
|---|---|---|---|
| Predictors | Odds Ratios | CI | p |
| (Intercept) | 1.28 | 0.67 – 2.42 | 0.456 |
| Treatment [Nuclei_25] | 0.55 | 0.23 – 1.33 | 0.185 |
| Treatment [Nuclei_50] | 0.51 | 0.21 – 1.24 | 0.136 |
| Treatment [Strip_25] | 1.23 | 0.52 – 2.91 | 0.645 |
| Treatment [Strip_50] | 1.28 | 0.54 – 3.02 | 0.572 |
| Planted [Outside] | 0.42 | 0.19 – 0.91 | 0.028 |
| Treatment [Nuclei_25] × Planted [Outside] |
1.14 | 0.37 – 3.52 | 0.824 |
| Treatment [Nuclei_50] × Planted [Outside] |
1.43 | 0.45 – 4.57 | 0.543 |
| Treatment [Strip_25] × Planted [Outside] |
0.46 | 0.15 – 1.42 | 0.178 |
| Random Effects | |||
| σ2 | 0.37 | ||
| τ00 Block | 0.00 | ||
| N Block | 5 | ||
| Observations | 176 | ||
| Marginal R2 / Conditional R2 | 0.500 / NA | ||
(iii) Only considering Baccharis + Vernonanthura
Table S2 - Only ruderal
Estimating density of recruits (\(ind.ha^{-1}\)) and without ruderal (Baccharis and Vernonanthura)
show code
density_ruderal_n_m2 <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
# dplyr::select(Treatment, Block, Plot, Planted)|>unique()|> with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
dplyr::add_row(
scientificName = c(NA, NA),
`Seed dispersal` = c(NA, NA),
`Planting group` = c(NA, NA),
Treatment = c("Full_Planted", "Full_Planted"),
Block = c("B5", "B1"),
Plot = c("P2", "P3"),
Planted = c("Inside","Inside"),
n_sapling= c(0,0),
n_tree = c(0,0)
) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop")
dens_ruderal_ha <- density_ruderal_n_m2 |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
)
dens_ruderal_P_UP <- dens_ruderal_ha |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
result_ruderal <- dens_ruderal_P_UP |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Smaller recruits
- Only in planted areas
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Inside, na.rm=TRUE),1),
min = round(min(n_small_ha_Inside, na.rm=TRUE),0),
max = round(max(n_small_ha_Inside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 2545.1 (1580.0 - 4566.0) 4388.9 (2847.0 - 6858.0) 5145.8 (1… 4093.8 … 5409.7 …- Only in unplanted areas
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
min = round(min(n_small_ha_Outside, na.rm=TRUE), 0),
max = round(max(n_small_ha_Outside, na.rm=TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_small_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 2784.7 (1181.0 - 6858.0) 2638.9 (1302.0 - … 2246.5 … 2489.6 …- In total (planted + unplanted areas)
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_small_total_ha, na.rm=TRUE),1),
min = round(min(n_small_total_ha, na.rm=TRUE),0),
max = round(max(n_small_total_ha, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 2545.1 (1580.0 - 4566.0) 1592.9 (712.0 - 5143.0) 1946.2 (65… 1354.2 … 1974.8 …Larger recruits
All treatments and inside/outside areas showed similar values
- Only in planted areas
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Inside, na.rm=TRUE),1),
min = round(min(n_large_ha_Inside, na.rm=TRUE),0),
max = round(max(n_large_ha_Inside, na.rm=TRUE),0)
) |>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 13.9 (0.0 - 35.0) 17.4 (0.0 - 69.0) 3.5 (0.0 - 17.0) 3.5 (0.0 - 17.0) 3.5 (0.…- Only in unplanted areas
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_ha_Outside, na.rm=TRUE),1),
min = round(min(n_large_ha_Outside, na.rm=TRUE),0),
max = round(max(n_large_ha_Outside, na.rm=TRUE),0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)Warning: There were 2 warnings in `dplyr::summarise()`.
The first warning was:
ℹ In argument: `min = round(min(n_large_ha_Outside, na.rm = TRUE), 0)`.
ℹ In group 1: `Treatment = "Full_Planted"`.
Caused by warning in `min()`:
! nenhum argumento não faltante para min; retornando Inf
ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 NaN (Inf - -Inf) 31.2 (0.0 - 52.0) 24.3 (0.0 - 87.0) 48.6 (0.0 - 139… 17.4 (0…- In total (planted + unplanted areas)
show code
result_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean = round(mean(n_large_total_ha, na.rm = TRUE), 1),
min = round(min(n_large_total_ha, na.rm = TRUE), 0),
max = round(max(n_large_total_ha, na.rm = TRUE), 0)
)|>
dplyr::mutate(
larger = sprintf("%.1f (%.1f - %.1f)",
mean,
min,
max)) |>
dplyr::select(Treatment, larger) |>
tidyr::pivot_wider(
names_from = Treatment,
values_from = larger
)# A tibble: 1 × 5
Full_Planted Nuclei_25 Nuclei_50 Strip_25 Strip_50
<chr> <chr> <chr> <chr> <chr>
1 13.9 (0.0 - 35.0) 13.9 (0.0 - 39.0) 6.9 (0.0 - 43.0) 18.7 (0.0 - 104… 5.2 (0.…Analysis for their abundances
show code
ruderal_df <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes"))|>
# dplyr::select(Treatment, Block, Plot, Planted)|>unique()|> with(table(Block, Plot, by=Treatment)) # Full_Planted P2B5 & P3B1 (Inside)
dplyr::add_row(
scientificName = c(NA, NA),
`Seed dispersal` = c(NA, NA),
`Planting group` = c(NA, NA),
Treatment = c("Full_Planted", "Full_Planted"),
Block = c("B5", "B1"),
Plot = c("P2", "P3"),
Planted = c("Inside","Inside"),
n_sapling= c(0,0),
n_tree = c(0,0),
n_total= c(0,0)
)|>
dplyr::group_by(Treatment, Block, Plot, Planted)|>
dplyr::summarise(n = sum(n_total), .groups="drop")Plotting
show code
ruderal_df |>
ggplot(aes(x=Planted, y=n, color=Block))+
geom_point()+facet_wrap(~Treatment)+theme_bw()
Models
show code
modelo_rud1 <- glmmTMB::glmmTMB(n ~ Treatment*Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
modelo_rud2 <- glmmTMB::glmmTMB(n ~ Treatment + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)
modelo_rud3 <- glmmTMB::glmmTMB(n ~ Planted + (1 | Block),
ziformula = ~1,
# family = poisson(),
family = glmmTMB::nbinom2(),
# family = glmmTMB::truncated_nbinom2(),
data = ruderal_df)
modelo_rud_null <- glmmTMB::glmmTMB(n ~ 0 + (1 | Block),
ziformula = ~1,
# family = poisson(), # ZIP
family = glmmTMB::nbinom2(), # ZINB
# family = glmmTMB::truncated_nbinom2(), # hurdle
data = ruderal_df)
# Model (null) | (treat*planted) | (treat) | (planted)
# ZIP
AIC(modelo_rud_null, modelo_rud1, modelo_rud2, modelo_rud3) df AIC
modelo_rud_null 3 1725.166
modelo_rud1 12 1658.610
modelo_rud2 8 1707.148
modelo_rud3 5 1660.457show code
performance::check_overdispersion(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside# Overdispersion test
dispersion ratio = 0.965
p-value = 0.872No overdispersion detected.show code
performance::check_model(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutside`check_outliers()` does not yet support models of class `glmmTMB`.
Exploring residuals
show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotQQunif(modelo_rud1)dropping columns from rank-deficient conditional model: TreatmentStrip_50:PlantedOutsideshow code
DHARMa::plotResiduals(modelo_rud3)
DHARMa::plotQQunif(modelo_rud3)
show code
car::Anova(modelo_rud1)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: n
Chisq Df Pr(>Chisq)
Treatment 16.1644 4 0.002806 **
Planted 65.9720 1 4.574e-16 ***
Treatment:Planted 1.5434 3 0.672287
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Is there any difference between treatments considering inside/outside?
show code
emmeans::emmeans(modelo_rud1, ~ Treatment|Planted) |>
multcomp::cld(adjust = "bonferroni")Planted = Inside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Full_Planted 3.66 0.206 Inf 3.13 4.19 1
Strip_25 4.09 0.201 Inf 3.57 4.60 12
Nuclei_25 4.13 0.201 Inf 3.61 4.65 12
Nuclei_50 4.23 0.200 Inf 3.72 4.75 2
Strip_50 4.29 0.200 Inf 3.78 4.81 2
Planted = Outside:
Treatment emmean SE df asymp.LCL asymp.UCL .group
Strip_25 3.42 0.203 Inf 2.91 3.93 1
Strip_50 3.48 0.202 Inf 2.97 3.99 1
Nuclei_50 3.51 0.202 Inf 3.00 4.01 1
Nuclei_25 3.60 0.204 Inf 3.09 4.11 1
Full_Planted nonEst NA NA NA NA
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
P value adjustment: bonferroni method for varying numbers of tests
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: Yes (planted) and no (only unplanted areas)
Is there any difference between inside/outside considering treatments?
show code
emmeans::emmeans(modelo_rud1, ~ Planted|Treatment) |>
multcomp::cld(adjust = "bonferroni")Treatment = Full_Planted:
Planted emmean SE df asymp.LCL asymp.UCL .group
Inside 3.66 0.206 Inf 3.26 4.06 1
Outside nonEst NA NA NA NA
Treatment = Nuclei_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.60 0.204 Inf 3.15 4.06 1
Inside 4.13 0.201 Inf 3.68 4.58 2
Treatment = Nuclei_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.51 0.202 Inf 3.05 3.96 1
Inside 4.23 0.200 Inf 3.79 4.68 2
Treatment = Strip_25:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.42 0.203 Inf 2.97 3.88 1
Inside 4.09 0.201 Inf 3.64 4.54 2
Treatment = Strip_50:
Planted emmean SE df asymp.LCL asymp.UCL .group
Outside 3.48 0.202 Inf 3.03 3.93 1
Inside 4.29 0.200 Inf 3.85 4.74 2
Results are given on the log (not the response) scale.
Confidence level used: 0.95
Conf-level adjustment: bonferroni method for varying numbers of estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Answer: No.
Table S3
show code
sjPlot::tab_model(modelo_rud1)Model matrix is rank deficient. Parameters
`TreatmentStrip_50:PlantedOutside` were not estimable.| n | ||||
|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | CI | p | |
| Count Model | ||||
| (Intercept) | 38.84 | 25.95 – 58.15 | <0.001 | |
| Treatment [Nuclei_25] | 1.60 | 1.14 – 2.24 | 0.006 | |
| Treatment [Nuclei_50] | 1.78 | 1.27 – 2.49 | 0.001 | |
| Treatment [Strip_25] | 1.53 | 1.09 – 2.15 | 0.013 | |
| Treatment [Strip_50] | 1.89 | 1.35 – 2.64 | <0.001 | |
| Planted [Outside] | 0.44 | 0.32 – 0.61 | <0.001 | |
| Treatment [Nuclei_25] × Planted [Outside] |
1.33 | 0.84 – 2.12 | 0.228 | |
| Treatment [Nuclei_50] × Planted [Outside] |
1.09 | 0.69 – 1.73 | 0.711 | |
| Treatment [Strip_25] × Planted [Outside] |
1.16 | 0.73 – 1.85 | 0.529 | |
| (Intercept) | 51.38 | 23.07 – 140.38 | ||
| Zero-Inflated Model | ||||
| (Intercept) | 0.01 | 0.00 – 0.05 | <0.001 | |
| Random Effects | ||||
| σ2 | 0.02 | |||
| τ00 Block | 0.13 | |||
| ICC | 0.87 | |||
| N Block | 5 | |||
| Observations | 179 | |||
| Marginal R2 / Conditional R2 | 0.424 / 0.924 | |||
show code
broom.mixed::tidy(modelo_rud1) |>
dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3))) |>
dplyr::select(-c(component, group))|>
knitr::kable()| effect | term | estimate | std.error | statistic | p.value |
|---|---|---|---|---|---|
| fixed | (Intercept) | 3.659 | 0.206 | 17.777 | 0.000 |
| fixed | TreatmentNuclei_25 | 0.471 | 0.171 | 2.748 | 0.006 |
| fixed | TreatmentNuclei_50 | 0.575 | 0.171 | 3.356 | 0.001 |
| fixed | TreatmentStrip_25 | 0.427 | 0.172 | 2.486 | 0.013 |
| fixed | TreatmentStrip_50 | 0.635 | 0.171 | 3.713 | 0.000 |
| fixed | PlantedOutside | -0.814 | 0.167 | -4.873 | 0.000 |
| fixed | TreatmentNuclei_25:PlantedOutside | 0.286 | 0.237 | 1.207 | 0.228 |
| fixed | TreatmentNuclei_50:PlantedOutside | 0.087 | 0.236 | 0.371 | 0.711 |
| fixed | TreatmentStrip_25:PlantedOutside | 0.149 | 0.237 | 0.630 | 0.529 |
| fixed | TreatmentStrip_50:PlantedOutside | NA | NA | NA | NA |
| fixed | (Intercept) | -4.497 | 0.721 | -6.238 | 0.000 |
| ran_pars | sd__(Intercept) | 0.365 | NA | NA | NA |
(iv) Figure Recruitment
Small letters: Comparison between treatments, only inside planted areas.
Capital letters: Comparison between treatments, only outside planted areas.
Smaller recruits
show code
ordem <- c("Nuclei_25", "Strip_25", "Nuclei_50", "Strip_50", "Full_Planted")
recruits_ha <- species_table |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop") |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
) |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
) |>
dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)
exc_recruits_ha <- species_table |>
dplyr::filter(!scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop") |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
) |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
) |>
dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)
ruderal_ha <- species_table |>
dplyr::filter(scientificName %in% c("Baccharis dracunculifolia", "Vernonanthura polyanthes")) |>
dplyr::mutate(
n_small = dplyr::if_else(is.na(n_sapling), 0, n_sapling),
n_large = dplyr::if_else(is.na(n_tree), 0, n_tree)) |>
dplyr::group_by(Treatment, Block, Planted) |>
dplyr::summarise(
n_small_m2 = (sum(n_small, na.rm=T))/(144*4),
n_large_m2 = (sum(n_large, na.rm=T))/(144*4),
.groups = "drop") |>
dplyr::mutate(
n_small_ha = n_small_m2 * 10000,
n_large_ha = n_large_m2 * 10000
) |>
tidyr::pivot_wider(
names_from = Planted,
values_from = c(n_small_ha, n_large_ha),
names_sep = "_"
) |>
dplyr::select(Treatment, Block, n_small_ha_Inside:n_large_ha_Outside)
prop_table <- tidyr::tibble(
Treatment = unique(sort(species_table$Treatment)),
prop_P = c(1, 0.25, 0.5, 0.25, 0.5),
prop_UP = 1 - prop_P
)
data_fig_recruits <- recruits_ha |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)
data_fig_ruderal <- ruderal_ha |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)
data_fig_exc_ruderal <- exc_recruits_ha |>
dplyr::left_join(prop_table, by = "Treatment") |>
dplyr::mutate(
n_small_total_ha = dplyr::coalesce(n_small_ha_Inside * prop_P, 0) + dplyr::coalesce(n_small_ha_Outside * prop_UP, 0),
n_large_total_ha = dplyr::coalesce(n_large_ha_Inside * prop_P, 0) + dplyr::coalesce(n_large_ha_Outside * prop_UP, 0)
) |>
dplyr::select(Treatment,
starts_with("n_small_ha"),
starts_with("n_large_ha"),
n_small_total_ha,
n_large_total_ha
)Small regenerants
show code
data_fig_small_1 <- data_fig_recruits |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(class="all recruits")
data_fig_small_2 <- data_fig_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(class="ruderal")
data_fig_small_3 <- data_fig_exc_ruderal |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_small_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_small_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_small_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_small_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::filter(!is.na(se)) |>
dplyr::mutate(class = "except ruderal")
data_fig_small_alterado <- data_fig_small_1
data_fig_small_alterado$mean <- data_fig_small_3$mean
data_fig_small <- dplyr::bind_rows(data_fig_small_alterado, data_fig_small_2)New Figure
show code
ordem_x <- c(
"Nuclei_25_inside", "Nuclei_25_outside",
"Strip_25_inside", "Strip_25_outside",
"Nuclei_50_inside", "Nuclei_50_outside",
"Strip_50_inside", "Strip_50_outside",
"Full_Planted_inside"
)
data_fig_small_top2 <- data_fig_small |>
group_by(Treatment, side) |>
arrange(desc(class)) |>
mutate(cum_mean = cumsum(mean)) |>
ungroup() |>
mutate(Treat_side = interaction(Treatment, side, sep = "_"),
Treat_side = factor(Treat_side, levels = ordem_x))|>
dplyr::arrange(desc(class),Treat_side)
letras <- data.frame(
Treat_side = levels(data_fig_small_top2$Treat_side),
label = c("AB","a",
"AB","a",
"B","a",
"B","a",
"A"),
y = c(5500, 5000, 5200, 4500,
7000, 4500, 7000, 4800, 3800)
)
fig_small<-ggplot(data_fig_small_top2,
aes(x = Treat_side, y = mean, fill = side, pattern = class)) +
geom_col_pattern(
position = "stack",
color = "black",
pattern_fill = "#92D050", pattern_alpha = 0.5,
pattern_angle = 45, pattern_density = 0.1,
pattern_spacing = 0.1
) +
scale_pattern_manual(
values = c("regenerant" = "none", "ruderal" = "stripe"),
guide = "none") +
scale_fill_manual(labels = c("Planted", "Unplanted"),
values = c("#92D050", "white")) +
geom_errorbar(
data = data_fig_small_top2 |> group_by(Treatment, side) |> slice_tail(n = 1),
aes(ymin = cum_mean - se, ymax = cum_mean + se, x = Treat_side),
width = .2,
inherit.aes = FALSE
) +
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}~")"),
x = "") +
geom_text(data = letras,
aes(x = Treat_side, y = y, label = label),
inherit.aes = FALSE,
size = 3.5, vjust = -0.5) +
scale_x_discrete(labels = function(x) gsub("_", "\n", x)) +
scale_y_continuous(expand =c(0,0), limits = c(0,8000), breaks = seq(0, 8000, by=1500))+
theme_classic() +
theme(
axis.text.x = element_blank(), #element_text(color="black", hjust = 0.5, size=9),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
legend.position = c(0.18, 0.95),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))
fig_small
show code
# ggsave("fig_small_recruits.png", plot=fig_small, width=90, height = 90 ,unit="mm")Larger recruits
show code
data_fig_large <- result_regenerant_all |>
dplyr::group_by(Treatment)|>
dplyr::summarise(
mean_inside = round(mean(n_large_ha_Inside, na.rm=TRUE), 1),
se_inside = round((sd(n_large_ha_Inside, na.rm=TRUE))/sqrt(dplyr::n()), 1),
mean_outside = round(mean(n_large_ha_Outside, na.rm=TRUE), 1),
se_outside = round((sd(n_large_ha_Outside, na.rm=TRUE))/sqrt(dplyr::n()), 1)
) |>
tidyr::pivot_longer(cols = -Treatment, names_to = c(".value", "side"), names_sep = "_") |>
dplyr::mutate(
Treatment = ordered(Treatment,levels = ordem),
tratamento = forcats::fct_recode(Treatment,
"Nuclei 25%" = "Nuclei_25", "Strips 25%" = "Strip_25",
"Nuclei 50%" = "Nuclei_50", "Strips 50%" = "Strip_50","Planting 100%" = "Full_Planted"),
Treat_side = interaction(Treatment, side, sep = "_"),
Treat_side = factor(Treat_side, ordered = TRUE, levels = ordem_x,
labels = c("Nuclei\n 25%", " " ,
"Strips\n 25%", " ",
"Nuclei\n 50%", " ",
"Strips\n 50%", " ",
"Planting\n100%"))
)|>
dplyr::filter(!is.na(se))Another way
show code
letras <- data.frame(
Treat_side = levels(data_fig_large$Treat_side),
label = c("A","a",
"A","a",
"A","a",
"A","a",
"A"),
y = c(310, 470,
200, 370,
200, 200,
300, 470,
320)
)
data_fig_large <- data_fig_large |>
dplyr::mutate(hjust_value = ifelse(Treat_side == last(Treat_side), 0.5, 0))
fig_large<-ggplot(data_fig_large,
aes(x = Treat_side, y = mean, fill = side)) +
geom_bar(
stat = "identity",
position = "stack",
color = "black") +
scale_fill_manual(labels = c("Planted", "Unplanted"),
values = c("#92D050", "white")) +
geom_errorbar(
aes(ymin = mean - se, ymax = mean + se, x = Treat_side),
width = .2
) +
labs(fill = "",
y = expression("Recruits (n.ha"^{-1}*")"),
x = "") +
geom_text(data = letras, aes(x = Treat_side, y = y, label = label), inherit.aes = FALSE, size = 3.5, vjust = -0.5) +
scale_y_continuous(expand =c(0,0), limits = c(0,550), breaks = seq(0, 500, by=100))+
theme_classic() +
theme(
axis.text.x = element_text(color="black", hjust = data_fig_large$hjust_value, size=9),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=11),
legend.position = "none",# c(0.18, 0.95),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
legend.spacing = unit(0, "mm"),
legend.background = element_blank(),
legend.title = element_text(size = 8))Warning: Vectorized input to `element_text()` is not officially supported.
ℹ Results may be unexpected or may change in future versions of ggplot2.show code
fig_large
Altogether:
show code
fig_recruits <- cowplot::plot_grid(fig_small + theme(axis.text.x=element_blank()),
fig_large,
ncol = 1,
labels = c("(a)", "(b)"), label_size = 10, label_x = 0.01, label_y = c(1,1.1)
)
fig_recruits
show code
ggsave("fig_recruits.png", plot=fig_recruits, width=90, height = 100 ,unit="mm")(vi) Species list table
Only recruits (not including planted individuals)
show code
table_recruits <- readxl::read_xlsx("Itatinga_Nuc_recruitment_SpeciesList.xlsx")
table_recruits |> knitr::kable(format = "html")| family | scientificName | Nuclei 25% (small) | Nuclei 25% (large) | Nuclei 25% (total) | Strips 25% (small) | Strips 25% (large) | Strips 25% (total) | Nuclei 50% (small) | Nuclei 50% (large) | Nuclei 50% (total) | Strips 50% (small) | Strips 50% (large) | Strips 50% (total) | Planting 100% (small) | Planting 100% (large) | Planting 100% (total) | Total small | Total large | Total |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| Anacardiaceae | Astronium graveolens | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 5 | 5 | 0 | 0 | 0 | 1 | 5 | 6 |
| Anacardiaceae | Myracrodruon urundeuva | 0 | 12 | 12 | 0 | 2 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 16 | 16 |
| Anacardiaceae | Schinus terebinthifolia | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Anacardiaceae | Tapirira guianensis | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Annonaceae | Annona cacans | 0 | 0 | 0 | 11 | 0 | 11 | 0 | 0 | 0 | 0 | 2 | 2 | 1 | 0 | 1 | 12 | 2 | 14 |
| Annonaceae | Annona coriacea | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 4 |
| Annonaceae | Duguetia furfuracea | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 1 | 2 | 1 | 3 |
| Apocynaceae | Aspidosperma tomentosum | 2 | 1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 | 2 | 5 | 7 |
| Araliaceae | Didymopanax vinosus | 0 | 0 | 0 | 4 | 3 | 7 | 2 | 0 | 2 | 3 | 1 | 4 | 0 | 0 | 0 | 9 | 4 | 13 |
| Arecaceae | Syagrus romanzoffiana | 1 | 1 | 2 | 1 | 0 | 1 | 3 | 0 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 1 | 6 |
| Asteraceae | Baccharis dracunculifolia | 1432 | 14 | 1446 | 916 | 15 | 931 | 931 | 8 | 939 | 1784 | 6 | 1790 | 511 | 4 | 515 | 5574 | 47 | 5621 |
| Asteraceae | Moquiniastrum barrosoae | 0 | 0 | 0 | 3 | 1 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 1 | 4 |
| Asteraceae | Moquiniastrum polymorphum | 0 | 17 | 17 | 2 | 3 | 5 | 13 | 2 | 15 | 2 | 0 | 2 | 0 | 0 | 0 | 17 | 22 | 39 |
| Asteraceae | Piptocarpha rotundifolia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Asteraceae | Vernonanthura polyanthes | 634 | 0 | 634 | 910 | 0 | 910 | 1311 | 0 | 1311 | 491 | 0 | 491 | 0 | 2 | 2 | 3346 | 2 | 3348 |
| Bignoniaceae | Cybistax antisyphilitica | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Bignoniaceae | Handroanthus albus | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 9 | 0 | 9 | 9 | 2 | 11 |
| Bignoniaceae | Handroanthus ochraceus | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 3 |
| Boraginaceae | Cordia sellowiana | 0 | 0 | 0 | 0 | 6 | 6 | 0 | 3 | 3 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 12 | 12 |
| Burseraceae | Protium spruceanum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Calophyllaceae | Myrsine gardneriana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Cannabaceae | Trema micrantha | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 10 | 10 | 0 | 0 | 0 | 0 | 11 | 11 |
| Chrysobalanaceae | Couepia grandiflora | 1 | 0 | 1 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 3 |
| Combretaceae | Terminalia argentea | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 |
| Euphorbiaceae | Alchornea triplinervia | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 44 | 0 | 44 | 45 | 0 | 45 |
| Euphorbiaceae | Croton floribundus | 73 | 0 | 73 | 55 | 1 | 56 | 119 | 0 | 119 | 147 | 0 | 147 | 0 | 0 | 0 | 394 | 1 | 395 |
| Euphorbiaceae | Croton urucurana | 1 | 0 | 1 | 2 | 0 | 2 | 0 | 0 | 0 | 7 | 1 | 8 | 0 | 0 | 0 | 10 | 1 | 11 |
| Euphorbiaceae | Mabea fistulifera | 1 | 0 | 1 | 12 | 0 | 12 | 7 | 0 | 7 | 31 | 2 | 33 | 0 | 3 | 3 | 51 | 5 | 56 |
| Euphorbiaceae | Sapium glandulosum | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 1 | 2 | 3 |
| Fabaceae | Anadenanthera colubrina var. cebil | 2 | 5 | 7 | 5 | 0 | 5 | 0 | 1 | 1 | 8 | 13 | 21 | 0 | 0 | 0 | 15 | 19 | 34 |
| Fabaceae | Andira fraxinifolia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Bauhinia rufa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 2 | 0 | 2 |
| Fabaceae | Chamaecrista comosa | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| Fabaceae | Copaifera langsdorffii | 2 | 2 | 4 | 2 | 4 | 6 | 2 | 2 | 4 | 5 | 0 | 5 | 0 | 0 | 0 | 11 | 8 | 19 |
| Fabaceae | Dalbergia frutescens | 0 | 2 | 2 | 7 | 0 | 7 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 7 | 2 | 9 |
| Fabaceae | Dalbergia miscolobium | 0 | 0 | 0 | 1 | 2 | 3 | 1 | 1 | 2 | 2 | 1 | 3 | 0 | 0 | 0 | 4 | 4 | 8 |
| Fabaceae | Enterolobium contortisiliquum | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Erythrina mulungu | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 |
| Fabaceae | Hymenaea stigonocarpa | 0 | 2 | 2 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Inga vera | 0 | 0 | 0 | 4 | 0 | 4 | 4 | 0 | 4 | 3 | 0 | 3 | 1 | 0 | 1 | 12 | 0 | 12 |
| Fabaceae | Leptolobium dasycarpum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Leptolobium elegans | 2 | 32 | 34 | 1 | 17 | 18 | 3 | 32 | 35 | 3 | 0 | 3 | 1 | 0 | 1 | 10 | 81 | 91 |
| Fabaceae | Leucochloron incuriale | 0 | 0 | 0 | 0 | 5 | 5 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 7 | 7 |
| Fabaceae | Machaerium acutifolium | 1 | 0 | 1 | 1 | 3 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 3 | 5 |
| Fabaceae | Machaerium brasiliense | 1 | 0 | 1 | 1 | 5 | 6 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 3 | 5 | 8 |
| Fabaceae | Machaerium villosum | 0 | 0 | 0 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 |
| Fabaceae | Platypodium elegans | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 2 | 0 | 2 |
| Fabaceae | Pterogyne nitens | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 |
| Fabaceae | Senegalia polyphylla | 48 | 0 | 48 | 94 | 1 | 95 | 87 | 2 | 89 | 101 | 0 | 101 | 0 | 0 | 0 | 330 | 3 | 333 |
| Fabaceae | Senna pendula | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Fabaceae | Senna rugosa | 6 | 0 | 6 | 7 | 0 | 7 | 9 | 0 | 9 | 7 | 35 | 42 | 0 | 28 | 28 | 29 | 63 | 92 |
| Lamiaceae | Aegiphila verticillata | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Lamiaceae | Vitex polygama | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 14 | 0 | 14 | 14 | 2 | 16 |
| Lauraceae | Aiouea sellowiana | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 4 |
| Lauraceae | Endlicheria paniculata | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 2 | 4 | 0 | 0 | 0 | 2 | 2 | 4 |
| Lauraceae | Ocotea corymbosa | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 3 |
| Lauraceae | Ocotea puberula | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Lauraceae | Ocotea pulchella | 9 | 5 | 14 | 2 | 0 | 2 | 2 | 1 | 3 | 0 | 6 | 6 | 0 | 0 | 0 | 13 | 12 | 25 |
| Lauraceae | Ocotea velutina | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 1 | 1 |
| Lythraceae | Lafoensia pacari | 0 | 0 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 |
| Malpighiaceae | Byrsonima crassifolia | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Malpighiaceae | Byrsonima intermedia | 0 | 0 | 0 | 6 | 1 | 7 | 3 | 0 | 3 | 0 | 6 | 6 | 0 | 0 | 0 | 9 | 7 | 16 |
| Malvaceae | Eriotheca gracilipes | 1 | 1 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 |
| Melastomataceae | Miconia albicans | 13 | 12 | 25 | 2 | 6 | 8 | 0 | 3 | 3 | 9 | 4 | 13 | 7 | 0 | 7 | 31 | 25 | 56 |
| Melastomataceae | Miconia auricoma | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 11 | 11 | 1 | 3 | 4 | 2 | 14 | 16 |
| Melastomataceae | Miconia langsdorffii | 0 | 5 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 2 | 5 | 7 |
| Melastomataceae | Miconia ligustroides | 5 | 0 | 5 | 5 | 0 | 5 | 0 | 0 | 0 | 6 | 7 | 13 | 0 | 2 | 2 | 16 | 9 | 25 |
| Melastomataceae | Miconia minutiflora | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 2 | 0 | 2 | 2 | 3 | 5 |
| Myrtaceae | Campomanesia pubescens | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 2 | 0 | 0 | 0 | 1 | 1 | 2 |
| Myrtaceae | Eugenia aurata | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 6 | 7 | 0 | 0 | 0 | 2 | 6 | 8 |
| Myrtaceae | Eugenia bimarginata | 1 | 0 | 1 | 1 | 2 | 3 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 3 | 2 | 5 |
| Myrtaceae | Eugenia hyemalis | 1 | 0 | 1 | 9 | 0 | 9 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 0 | 10 |
| Myrtaceae | Eugenia myrcianthes | 0 | 0 | 0 | 2 | 1 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 1 | 3 |
| Myrtaceae | Eugenia punicifolia | 5 | 0 | 5 | 6 | 0 | 6 | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 12 | 0 | 12 |
| Myrtaceae | Eugenia pyriformis | 1 | 0 | 1 | 6 | 0 | 6 | 0 | 0 | 0 | 6 | 0 | 6 | 0 | 0 | 0 | 13 | 0 | 13 |
| Myrtaceae | Eugenia uniflora | 0 | 0 | 0 | 5 | 0 | 5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 5 |
| Myrtaceae | Myrcia | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Myrtaceae | Myrcia bella | 2 | 14 | 16 | 3 | 11 | 14 | 5 | 8 | 13 | 2 | 3 | 5 | 0 | 0 | 0 | 12 | 36 | 48 |
| Myrtaceae | Myrcia guianensis | 5 | 5 | 10 | 1 | 0 | 1 | 4 | 2 | 6 | 0 | 0 | 0 | 0 | 0 | 0 | 10 | 7 | 17 |
| Myrtaceae | Myrcia rufipes | 1 | 4 | 5 | 2 | 0 | 2 | 1 | 0 | 1 | 3 | 0 | 3 | 0 | 0 | 0 | 7 | 4 | 11 |
| Myrtaceae | Myrcia splendens | 4 | 0 | 4 | 4 | 0 | 4 | 0 | 1 | 1 | 4 | 3 | 7 | 0 | 0 | 0 | 12 | 4 | 16 |
| Myrtaceae | Myrcia tomentosa | 0 | 0 | 0 | 0 | 2 | 2 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 3 | 3 |
| Myrtaceae | Myrciaria floribunda | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 6 | 6 | 0 | 0 | 0 | 1 | 6 | 7 |
| Myrtaceae | Myrciaria glazioviana | 11 | 0 | 11 | 4 | 1 | 5 | 5 | 0 | 5 | 7 | 0 | 7 | 0 | 0 | 0 | 27 | 1 | 28 |
| Myrtaceae | Psidium grandifolium | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 1 | 2 | 0 | 2 |
| Myrtaceae | Psidium guineense | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Myrtaceae | Psidium oligospermum | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 |
| Nyctaginaceae | Guapira noxia | 0 | 4 | 4 | 0 | 1 | 1 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 3 | 5 | 8 |
| Ochnaceae | Ouratea spectabilis | 9 | 12 | 21 | 4 | 9 | 13 | 1 | 0 | 1 | 6 | 0 | 6 | 0 | 0 | 0 | 20 | 21 | 41 |
| Peraceae | Pera glabrata | 0 | 0 | 0 | 1 | 0 | 1 | 4 | 8 | 12 | 1 | 1 | 2 | 0 | 2 | 2 | 6 | 11 | 17 |
| Piperaceae | Piper mollicomum | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Polygalaceae | Bredemeyera floribunda | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Polygonaceae | Coccoloba cordata | 0 | 0 | 0 | 3 | 5 | 8 | 1 | 0 | 1 | 1 | 0 | 1 | 0 | 0 | 0 | 5 | 5 | 10 |
| Polygonaceae | Triplaris americana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Primulaceae | Myrsine coriacea | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 6 | 8 | 0 | 0 | 0 | 2 | 6 | 8 |
| Primulaceae | Myrsine gardneriana | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 1 |
| Primulaceae | Myrsine guianensis | 0 | 0 | 0 | 4 | 0 | 4 | 0 | 3 | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 3 | 7 |
| Primulaceae | Myrsine umbellata | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 5 | 0 | 5 |
| Proteaceae | Roupala montana | 1 | 2 | 3 | 0 | 3 | 3 | 0 | 0 | 0 | 3 | 0 | 3 | 0 | 0 | 0 | 4 | 5 | 9 |
| Rubiaceae | Amaioua intermedia | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 4 | 0 | 4 | 18 | 0 | 18 | 23 | 0 | 23 |
| Rubiaceae | Palicourea sessilis | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 4 | 1 | 5 | 0 | 0 | 0 | 5 | 1 | 6 |
| Salicaceae | Casearia sylvestris | 36 | 0 | 36 | 76 | 0 | 76 | 60 | 0 | 60 | 15 | 0 | 15 | 2 | 0 | 2 | 189 | 0 | 189 |
| Sapindaceae | Matayba elaeagnoides | 0 | 0 | 0 | 3 | 0 | 3 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 5 | 0 | 5 |
| Sapindaceae | Sapindus saponaria | 0 | 0 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 3 | 3 | 11 | 0 | 11 | 12 | 3 | 15 |
| Sapotaceae | Pouteria ramiflora | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 1 | 1 | 26 | 0 | 26 | 26 | 2 | 28 |
| Solanaceae | Cestrum corymbosum | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 | 3 | 3 | 0 | 3 | 3 | 3 | 6 |
| Solanaceae | Solanum didymum | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 |
| Solanaceae | Solanum granuloso-leprosum | 44 | 3 | 47 | 79 | 6 | 85 | 77 | 1 | 78 | 55 | 3 | 58 | 2 | 7 | 9 | 257 | 20 | 277 |
| Solanaceae | Solanum sisymbriifolium | 166 | 0 | 166 | 72 | 0 | 72 | 80 | 0 | 80 | 98 | 0 | 98 | 0 | 0 | 0 | 416 | 0 | 416 |
| Urticaceae | Cecropia pachystachya | 1 | 0 | 1 | 2 | 0 | 2 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 9 | 9 | 4 | 9 | 13 |
| Verbenaceae | Lantana camara | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 0 | 12 | 12 | 1 | 0 | 1 | 3 | 12 | 15 |
| Vochysiaceae | Qualea grandiflora | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0 | 2 | 0 | 0 | 0 | 222 | 0 | 222 | 224 | 0 | 224 |
| NA | Unclassified | 9 | 8 | 17 | 5 | 4 | 9 | 8 | 1 | 9 | 12 | 2 | 14 | 0 | 0 | 0 | 34 | 15 | 49 |
! Altogether
show code
sampled_species <-readxl::read_xlsx("itatinga_sampled_species.xlsx")
sampled_species |>
dplyr::select(-`Seed dispersal`) |>
knitr::kable(format = "html")| family | scientificName | scientificNameAuthorship | Large_recruit | Small_recruit | Origin | Planting group | Successional group |
|---|---|---|---|---|---|---|---|
| Anacardiaceae | Astronium graveolens | Jacq. | TRUE | TRUE | planted | diversity | non-pioneer |
| Anacardiaceae | Schinus terebinthifolia | Raddi | FALSE | TRUE | recruit | NA | pioneer |
| Anacardiaceae | Myracrodruon urundeuva | Allem. | TRUE | FALSE | recruit | NA | non-pioneer |
| Anacardiaceae | Tapirira guianensis | Aubl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Annonaceae | Annona cacans | Warm. | TRUE | TRUE | planted | diversity | non-pioneer |
| Annonaceae | Annona coriacea | Mart. | TRUE | TRUE | recruit | NA | non-pioneer |
| Annonaceae | Duguetia furfuracea | (A.St.-Hil.) Saff. | TRUE | TRUE | recruit | NA | non-pioneer |
| Apocynaceae | Aspidosperma tomentosum | Mart. | TRUE | TRUE | recruit | NA | non-pioneer |
| Araliaceae | Didymopanax vinosus | (Cham. & Schltdl.) Marchal | TRUE | TRUE | recruit | NA | non-pioneer |
| Arecaceae | Syagrus romanzoffiana | (Cham.) Glassman | TRUE | TRUE | recruit | NA | non-pioneer |
| Asteraceae | Piptocarpha rotundifolia | Baker | FALSE | TRUE | recruit | NA | pioneer |
| Asteraceae | Baccharis dracunculifolia | DC. | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Moquiniastrum barrosoae | (Cabrera) G.Sancho | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Moquiniastrum polymorphum | (Less.) G.Sancho | TRUE | TRUE | recruit | NA | pioneer |
| Asteraceae | Vernonanthura polyanthes | (Spreng.) A.J.Vega & Dematt. | TRUE | TRUE | recruit | NA | pioneer |
| Bignoniaceae | Handroanthus ochraceus | (Cham.) Mattos | TRUE | TRUE | planted | diversity | non-pioneer |
| Bignoniaceae | Handroanthus impetiginosus | (Mart. ex DC.) Mattos | FALSE | FALSE | planted | diversity | non-pioneer |
| Bignoniaceae | Jacaranda cuspidifolia | Mart. | FALSE | FALSE | planted | diversity | non-pioneer |
| Bignoniaceae | Cybistax antisyphilitica | (Mart.) Mart. | TRUE | FALSE | recruit | NA | non-pioneer |
| Bignoniaceae | Handroanthus albus | (Cham.) Mattos | TRUE | TRUE | recruit | NA | non-pioneer |
| Boraginaceae | Cordia ecalyculata | Vell. | FALSE | FALSE | planted | diversity | pioneer |
| Boraginaceae | Cordia sellowiana | Cham. | TRUE | FALSE | recruit | NA | pioneer |
| Burseraceae | Protium spruceanum | Engl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Cannabaceae | Trema micrantha | (L.) Blume | TRUE | FALSE | recruit | NA | pioneer |
| Chrysobalanaceae | Couepia grandiflora | (Mart. & Zucc.) Benth. ex Hook.f. | TRUE | TRUE | recruit | NA | non-pioneer |
| Combretaceae | Terminalia argentea | Mart. | TRUE | FALSE | recruit | NA | non-pioneer |
| Euphorbiaceae | Croton floribundus | Spreng. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Croton urucurana | Baill. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Mabea fistulifera | Mart. | TRUE | TRUE | planted | filling | pioneer |
| Euphorbiaceae | Alchornea triplinervia | (Spreng.) Müll.Arg. | FALSE | TRUE | recruit | NA | non-pioneer |
| Euphorbiaceae | Sapium glandulosum | (L.) Morong | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Enterolobium contortisiliquum | (Vell.) Morong | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Erythrina mulungu | Mart. ex Benth. | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Platypodium elegans | Vogel | FALSE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Copaifera langsdorffii | Desf. | TRUE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Machaerium acutifolium | Vogel | TRUE | TRUE | planted | diversity | non-pioneer |
| Fabaceae | Pterogyne nitens | Tul. | TRUE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Anadenanthera colubrina var. cebil | (Griseb.) Altschul | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Anadenanthera peregrina | (L.) Speg. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Apuleia leiocarpa | (Vogel) J.F.Macbr. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Cassia fistula | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Hymenaea courbaril | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Myroxylon peruiferum | L.f. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Parapiptadenia rigida | (Benth.) Brenan | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Peltophorum dubium | (Spreng.) Taub. | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Senna macranthera | (DC. ex Collad.) H.S.Irwin & Barneby | FALSE | FALSE | planted | diversity | non-pioneer |
| Fabaceae | Inga vera | Willd. | FALSE | TRUE | planted | filling | non-pioneer |
| Fabaceae | Senegalia polyphylla | (DC.) Britton & Rose | TRUE | TRUE | planted | filling | pioneer |
| Fabaceae | Inga edulis | Mart. | FALSE | FALSE | planted | filling | non-pioneer |
| Fabaceae | Andira fraxinifolia | Benth. | FALSE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Bauhinia rufa | (Bong.) Steud. | FALSE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Chamaecrista comosa | E.Mey. | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Dalbergia frutescens | (Vell.) Britton | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Dalbergia miscolobium | Benth. | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Hymenaea stigonocarpa | Mart. ex Hayne | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Leptolobium dasycarpum | Vogel | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Leptolobium elegans | Vogel | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Leucochloron incuriale | (Vell.) Barneby & J.W.Grimes | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Machaerium brasiliense | Vogel | TRUE | TRUE | recruit | NA | non-pioneer |
| Fabaceae | Machaerium villosum | Vogel | TRUE | FALSE | recruit | NA | non-pioneer |
| Fabaceae | Senna pendula | (Humb. & Bonpl. ex Willd.) H.S.Irwin & Barneby | TRUE | FALSE | recruit | NA | pioneer |
| Fabaceae | Senna rugosa | (G.Don) H.S.Irwin & Barneby | TRUE | TRUE | recruit | NA | pioneer |
| Lamiaceae | Aegiphila verticillata | Vell. | FALSE | TRUE | recruit | NA | pioneer |
| Lamiaceae | Vitex polygama | Cham. | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea puberula | (Rich.) Nees | FALSE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Aiouea sellowiana | (Nees & Mart.) R.Rohde | TRUE | FALSE | recruit | NA | non-pioneer |
| Lauraceae | Endlicheria paniculata | (Spreng.) J.F.Macbr. | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea corymbosa | (Meisn.) Mez | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea pulchella | (Nees & Mart.) Mez | TRUE | TRUE | recruit | NA | non-pioneer |
| Lauraceae | Ocotea velutina | (Nees) Mart. ex B.D.Jacks. | TRUE | FALSE | recruit | NA | non-pioneer |
| Lecythidaceae | Cariniana estrellensis | (Raddi) Kuntze | FALSE | FALSE | planted | diversity | non-pioneer |
| Lythraceae | Lafoensia pacari | A.St.-Hil. | TRUE | TRUE | planted | diversity | non-pioneer |
| Malpighiaceae | Byrsonima crassifolia | Kunth | FALSE | TRUE | recruit | NA | non-pioneer |
| Malpighiaceae | Byrsonima intermedia | A.Juss. | TRUE | TRUE | recruit | NA | non-pioneer |
| Malvaceae | Ceiba speciosa | (A.St.-Hil., A.Juss. & Cambess.) Ravenna | FALSE | FALSE | planted | diversity | non-pioneer |
| Malvaceae | Heliocarpus popayanensis | Kunth | FALSE | FALSE | planted | diversity | pioneer |
| Malvaceae | Luehea divaricata | Mart. | FALSE | FALSE | planted | diversity | non-pioneer |
| Malvaceae | Guazuma ulmifolia | Lam. | FALSE | FALSE | planted | filling | pioneer |
| Malvaceae | Eriotheca gracilipes | (K.Schum.) A.Robyns | TRUE | TRUE | recruit | NA | non-pioneer |
| Melastomataceae | Miconia albicans | (Sw.) Steud. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia auricoma | (Spring ex Mart.) R.Goldenb. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia langsdorffii | Cogn. | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia ligustroides | Naudin | TRUE | TRUE | recruit | NA | pioneer |
| Melastomataceae | Miconia minutiflora | (Bonpl.) DC. | TRUE | TRUE | recruit | NA | pioneer |
| Meliaceae | Cedrela fissilis | Vell. | FALSE | FALSE | planted | diversity | non-pioneer |
| Moraceae | Ficus guaranitica | Chodat | FALSE | FALSE | planted | diversity | non-pioneer |
| Myrtaceae | Eugenia uniflora | L. | FALSE | TRUE | planted | diversity | non-pioneer |
| Myrtaceae | Psidium myrtoides | O.Berg | FALSE | FALSE | planted | diversity | non-pioneer |
| Myrtaceae | Eugenia hyemalis | Cambess. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia punicifolia | (Kunth) DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia pyriformis | Cambess. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia | DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium grandifolium | Mart. ex DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium guineense | Sw. | FALSE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Campomanesia pubescens | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia aurata | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia bimarginata | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Eugenia myrcianthes | Nied. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia bella | Cambess. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia guianensis | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia rufipes | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia splendens | DC. | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrcia tomentosa | DC. | TRUE | FALSE | recruit | NA | non-pioneer |
| Myrtaceae | Myrciaria floribunda | O.Berg | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Myrciaria glazioviana | (Kiaersk.) G.M.Barroso ex Sobral | TRUE | TRUE | recruit | NA | non-pioneer |
| Myrtaceae | Psidium oligospermum | Mart. ex DC. | TRUE | FALSE | recruit | NA | non-pioneer |
| Nyctaginaceae | Guapira noxia | (Netto) Lundell | TRUE | TRUE | recruit | NA | non-pioneer |
| Ochnaceae | Ouratea spectabilis | Engl. | TRUE | TRUE | recruit | NA | non-pioneer |
| Peraceae | Pera glabrata | (Schott) Poepp. ex Baill. | TRUE | TRUE | recruit | NA | non-pioneer |
| Petiveriaceae | Gallesia integrifolia | (Spreng.) Harms | FALSE | FALSE | planted | diversity | non-pioneer |
| Piperaceae | Piper mollicomum | Kunth | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygalaceae | Bredemeyera floribunda | Willd. | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygonaceae | Triplaris americana | L. | FALSE | TRUE | recruit | NA | non-pioneer |
| Polygonaceae | Coccoloba cordata | Cham. | TRUE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine gardneriana | A.DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine gardneriana | A.DC. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine umbellata | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine coriacea | (Sw.) R.Br. ex Roem. & Schult. | TRUE | TRUE | recruit | NA | non-pioneer |
| Primulaceae | Myrsine guianensis | (Aubl.) Kuntze | TRUE | TRUE | recruit | NA | non-pioneer |
| Proteaceae | Roupala montana | Aubl. | TRUE | TRUE | recruit | NA | non-pioneer |
| Rosaceae | Prunus myrtifolia | (L.) Urb. | FALSE | FALSE | planted | diversity | non-pioneer |
| Rubiaceae | Genipa americana | L. | FALSE | FALSE | planted | diversity | non-pioneer |
| Rubiaceae | Amaioua intermedia | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
| Rubiaceae | Palicourea sessilis | (Vell.) C.M.Taylor | TRUE | TRUE | recruit | NA | non-pioneer |
| Salicaceae | Casearia sylvestris | Sw. | FALSE | TRUE | recruit | NA | pioneer |
| Sapindaceae | Matayba elaeagnoides | Radlk. | FALSE | TRUE | recruit | NA | non-pioneer |
| Sapindaceae | Sapindus saponaria | L. | TRUE | TRUE | recruit | NA | non-pioneer |
| Sapotaceae | Pouteria ramiflora | Radlk. | TRUE | TRUE | recruit | NA | non-pioneer |
| Solanaceae | Solanum didymum | Dunal | FALSE | TRUE | recruit | NA | pioneer |
| Solanaceae | Solanum sisymbriifolium | Lam. | FALSE | TRUE | recruit | NA | pioneer |
| Solanaceae | Cestrum corymbosum | Schltdl. | TRUE | FALSE | recruit | NA | non-pioneer |
| Solanaceae | Solanum granuloso-leprosum | Dunal | TRUE | TRUE | recruit | NA | pioneer |
| Unclassified | Unclassified | NA | TRUE | TRUE | recruit | NA | NA |
| Urticaceae | Cecropia pachystachya | Trécul | TRUE | TRUE | planted | diversity | pioneer |
| Verbenaceae | Citharexylum myrianthum | Cham. | FALSE | FALSE | planted | diversity | non-pioneer |
| Verbenaceae | Lantana camara | L. | TRUE | TRUE | recruit | NA | pioneer |
| Vochysiaceae | Qualea grandiflora | Mart. | FALSE | TRUE | recruit | NA | non-pioneer |
D) Total AGB
(i) Organizing data
Estimation in Mg.ha-1
576 = 12m x 12m x 4 plots
show code
AGB_data0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, plantio)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |> # writexl::write_xlsx("AGB_plot.xlsx")
dplyr::mutate(
AGB_Mgha = dplyr::case_when(
tratamento == "100" ~ (AGB_sum / 576) * 10000,
tratamento %in% c("F50", "N50") & plantio == "CP" ~ (AGB_sum / 576) * 5000,
tratamento %in% c("F50", "N50") & plantio == "SP" ~ (AGB_sum / 576) * 5000,
tratamento %in% c("F25", "N25") & plantio == "CP" ~ (AGB_sum / 576) * 2500,
tratamento %in% c("F25", "N25") & plantio == "SP" ~ (AGB_sum / 576) * 7500),
Tratamento = ordered(tratamento, levels = ordem,
labels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")))
AGB_data <- AGB_data0 |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_Mgha = sum(AGB_Mgha), .groups = "drop") |>
dplyr::mutate(
Ano = ano-2021,
ANO = factor(ano),
scaled_AGB = scale(AGB_Mgha))(ii) Analysing Total aboveground biomass per hectare
Comparing models
show code
AGB_model0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_data)
AGB_model1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_data)
AGB_model2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_data) Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradientWarning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determinedshow code
AGB_model3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_data) # consider rescaling
AGB_model4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_data) # isSingular
AGB_model5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_data) # best model
AGB_model6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_data)
AGB_model7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_data) # isSingularboundary (singular) fit: see help('isSingular')show code
AIC(AGB_model0, AGB_model1, AGB_model2, AGB_model3, AGB_model4, AGB_model5,AGB_model6, AGB_model7 ) df AIC
AGB_model0 2 144.60044
AGB_model1 4 141.90535
AGB_model2 8 144.14482
AGB_model3 16 93.47468
AGB_model4 12 88.21817
AGB_model5 8 85.97871
AGB_model6 12 91.23521
AGB_model7 14 90.49850Exploring residuals
show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_model5)
DHARMa::plotQQunif(AGB_model5)
One outlier
Marginal \(R^2\) / Conditional \(R^2\)
show code
MuMIn::r.squaredGLMM(AGB_model5) # without interaction R2m R2c
[1,] 0.6955092 0.8374889show code
0.696/ 0.837[1] 0.8315412show code
MuMIn::r.squaredGLMM(AGB_model6) # with interaction R2m R2c
[1,] 0.701532 0.8419478show code
0.702/ 0.842[1] 0.83372920.702/ 0.842 - variance explained by fixed effects and by the entire model (around r0.702/ 0.842 % is explained by fixed effects)
show code
car::Anova(AGB_model6)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 165.4367 4 < 2.2e-16 ***
ANO 47.4491 1 5.645e-12 ***
tratamento:ANO 4.6061 4 0.3302
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction is no significant at all (\(p-value \sim 0.33\))
show code
car::Anova(AGB_model5)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 162.968 4 < 2.2e-16 ***
ANO 46.741 1 8.102e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Summary
show code
sjPlot::tab_model(AGB_model5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.96 | 0.51 – 1.42 | <0.001 |
| tratamento [F25] | -1.83 | -2.20 – -1.46 | <0.001 |
| tratamento [F50] | -1.35 | -1.72 – -0.98 | <0.001 |
| tratamento [N25] | -2.19 | -2.56 – -1.82 | <0.001 |
| tratamento [N50] | -1.44 | -1.82 – -1.07 | <0.001 |
| ANO [2024] | 0.80 | 0.56 – 1.03 | <0.001 |
| Random Effects | |||
| σ2 | 0.17 | ||
| τ00 bloco | 0.15 | ||
| ICC | 0.47 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.696 / 0.837 | ||
show code
broom.mixed::tidy(AGB_model5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 8 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.965 0.224 4.31
2 fixed <NA> tratamentoF25 -1.83 0.184 -9.92
3 fixed <NA> tratamentoF50 -1.35 0.184 -7.34
4 fixed <NA> tratamentoN25 -2.19 0.184 -11.9
5 fixed <NA> tratamentoN50 -1.44 0.184 -7.83
6 fixed <NA> ANO2024 0.797 0.117 6.84
7 ran_pars bloco sd__(Intercept) 0.385 NA NA
8 ran_pars Residual sd__Observation 0.412 NA NA Assessing comparisons:
show code
emmeans::emmeans(AGB_model6, pairwise ~ tratamento | ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
N25 -1.028 0.251 13.4 -1.5696 -0.4868 1
F25 -0.804 0.251 13.4 -1.3450 -0.2622 1
N50 -0.527 0.251 13.4 -1.0680 0.0148 1
F50 -0.423 0.251 13.4 -0.9645 0.1183 1
100 0.790 0.251 13.4 0.2484 1.3312 2
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.629 0.251 13.4 -1.1701 -0.0873 1
F25 -0.125 0.251 13.4 -0.6663 0.4165 12
N50 0.366 0.251 13.4 -0.1756 0.9072 2
F50 0.443 0.251 13.4 -0.0981 0.9847 2
100 1.936 0.251 13.4 1.3949 2.4777 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
AGB_model6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_data)
AIC(AGB_model6a, AGB_model6) df AIC
AGB_model6a 12 98.16669
AGB_model6 12 91.23521show code
emmeans::emtrends(AGB_model6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld() tratamento Ano.trend SE df lower.CL upper.CL .group
N25 0.200 0.129 36 -0.0625 0.462 1
F25 0.339 0.129 36 0.0771 0.602 1
F50 0.433 0.129 36 0.1709 0.695 1
N50 0.446 0.129 36 0.1840 0.708 1
100 0.573 0.129 36 0.3110 0.836 1
Results are averaged over the levels of: Ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. (iii) AGB only inside
Preparing data:
Comparison only inside plots over the years
show code
ordem <-c("N25", "F25", "N50", "F50", "100")
AGB_inside0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T),
.groups = "drop") |>
dplyr::mutate(
tratamento = ordered(tratamento, levels = ordem),
Parcela = paste(tratamento, bloco, parcela, sep="_"),
Ano = ano-2021) |>
dplyr::filter(plantio == "CP")View
show code
AGB_inside0 |>
ggplot(aes(x = factor(ano), y = log(AGB_sum), group = Parcela, colour = parcela)) +
geom_point() +
geom_line()+facet_grid(bloco~tratamento)+
theme_minimal()
Plots with decreasing :
show code
deltas <- AGB_inside0 |>
dplyr::select(ano, Parcela, AGB_sum)|>
tidyr::pivot_wider(names_from = ano, names_prefix = "ano_", values_from = AGB_sum)|>
dplyr::mutate(delta_AGB = (ano_2024/ano_2022))|>
dplyr::arrange(delta_AGB)
deltas |> head(10)# A tibble: 10 × 4
Parcela ano_2022 ano_2024 delta_AGB
<chr> <dbl> <dbl> <dbl>
1 100_4_4 422. 216. 0.511
2 F50_5_3 184. 128. 0.698
3 N50_5_4 324. 228. 0.701
4 F25_5_1 269. 198. 0.736
5 F50_1_1 152. 141. 0.930
6 N25_3_4 165. 157. 0.954
7 100_5_3 230. 220. 0.957
8 F50_3_4 208. 208. 0.998
9 N50_1_4 485. 493. 1.02
10 F50_2_4 549. 577. 1.05 show code
deltas |> tail(10)# A tibble: 10 × 4
Parcela ano_2022 ano_2024 delta_AGB
<chr> <dbl> <dbl> <dbl>
1 F50_2_1 216. 758. 3.51
2 N50_5_1 170. 611. 3.58
3 N50_5_3 153. 554. 3.62
4 N25_3_2 223. 846. 3.79
5 F50_2_3 77.5 313. 4.03
6 F50_2_2 82.8 367. 4.43
7 F50_3_3 175. 799. 4.56
8 F25_3_3 114. 557. 4.88
9 N50_4_2 136. 784. 5.75
10 N50_2_4 88.6 690. 7.79How much % is planted or recruit in AGB?
show code
trees_clean |>
dplyr::filter(plantio == "CP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T)/1000,
AGB_planted = sum(AGB_ferez[RP == "planted"], na.rm = TRUE)/1000,
prop_AGB_planted = (AGB_planted / AGB_sum)*100,
.groups = "drop") |> #dplyr::group_by(tratamento,plantio)|>
dplyr::summarise(mean(prop_AGB_planted),
sd(prop_AGB_planted)/sqrt(dplyr::n()))# A tibble: 1 × 2
`mean(prop_AGB_planted)` `sd(prop_AGB_planted)/sqrt(dplyr::n())`
<dbl> <dbl>
1 96.9 0.758Should I treat them with a random factor (1|bloco)?
show code
AGB_inside <- trees_clean |>
dplyr::filter(plantio == "CP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
dplyr::mutate(
Tratamento = ordered(tratamento, levels = ordem),
Ano = ano-2021,
ANO = ordered(ano),
scaled_AGB = scale(AGB_sum)) Comparing models
show code
AGB_inside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_inside)
AGB_inside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_inside)
AGB_inside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_inside) # best model Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradientWarning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determinedshow code
AGB_inside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_inside)
AGB_inside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_inside)
AGB_inside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_inside)
AGB_inside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_inside)
AGB_inside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_inside) # isSingularboundary (singular) fit: see help('isSingular')show code
AIC(AGB_inside_m0, AGB_inside_m1, AGB_inside_m2, AGB_inside_m3, AGB_inside_m4,AGB_inside_m5, AGB_inside_m6, AGB_inside_m7) df AIC
AGB_inside_m0 2 138.7732
AGB_inside_m1 4 117.1008
AGB_inside_m2 8 116.0631
AGB_inside_m3 16 122.2545
AGB_inside_m4 12 115.0812
AGB_inside_m5 8 116.1189
AGB_inside_m6 12 123.2922
AGB_inside_m7 14 123.2583Exploring residuals
show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_inside_m6)
DHARMa::plotQQunif(AGB_inside_m6)
Marginal \(R^2\) / Conditional \(R^2\)
show code
MuMIn::r.squaredGLMM(AGB_inside_m6) R2m R2c
[1,] 0.396267 0.699652show code
0.396/ 0.6996[1] 0.56603770.396/ 0.6996 variance explained by fixed effects and by the entire model (around 0.5660377 % is explained by fixed effects)
show code
car::Anova(AGB_inside_m6)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 14.3671 4 0.006211 **
ANO 47.6142 1 5.189e-12 ***
tratamento:ANO 2.6673 4 0.614946
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction is no significant at all (\(p-value \sim 0.85\))
show code
sjPlot::tab_model(AGB_inside_m5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.02 | -0.62 – 0.66 | 0.949 |
| tratamento [F25] | 0.50 | -0.02 – 1.01 | 0.059 |
| tratamento [F50] | 0.00 | -0.51 – 0.52 | 0.987 |
| tratamento [N25] | -0.47 | -0.99 – 0.04 | 0.070 |
| tratamento [N50] | -0.13 | -0.64 – 0.39 | 0.622 |
| ANO [linear] | 0.80 | 0.57 – 1.03 | <0.001 |
| Random Effects | |||
| σ2 | 0.33 | ||
| τ00 bloco | 0.34 | ||
| ICC | 0.51 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.390 / 0.702 | ||
show code
broom.mixed::tidy(AGB_inside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 8 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.02 0.318 0.064
2 fixed <NA> tratamentoF25 0.496 0.256 1.94
3 fixed <NA> tratamentoF50 0.004 0.256 0.016
4 fixed <NA> tratamentoN25 -0.475 0.256 -1.86
5 fixed <NA> tratamentoN50 -0.127 0.256 -0.497
6 fixed <NA> ANO.L 0.802 0.114 7.02
7 ran_pars bloco sd__(Intercept) 0.585 NA NA
8 ran_pars Residual sd__Observation 0.571 NA NA Assessing comparisons:
show code
emmeans::emmeans(AGB_inside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.8983 0.369 12.2 -1.700 -0.0969 1
N50 -0.7449 0.369 12.2 -1.546 0.0564 1
F50 -0.5376 0.369 12.2 -1.339 0.2637 1
100 -0.3908 0.369 12.2 -1.192 0.4105 1
F25 -0.2642 0.369 12.2 -1.066 0.5371 1
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N25 -0.0107 0.369 12.2 -0.812 0.7906 1
100 0.4316 0.369 12.2 -0.370 1.2329 12
N50 0.5318 0.369 12.2 -0.269 1.3332 12
F50 0.5866 0.369 12.2 -0.215 1.3879 12
F25 1.2965 0.369 12.2 0.495 2.0979 2
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
AGB_inside_m6a <- lme4::lmer(scaled_AGB ~ tratamento*Ano + (1|bloco), data = AGB_inside)
AIC(AGB_inside_m6a, AGB_inside_m6) df AIC
AGB_inside_m6a 12 126.7579
AGB_inside_m6 12 123.2922show code
emmeans::emtrends(AGB_inside_m6a, pairwise ~ tratamento, var = "Ano", poly.trend = TRUE) |> multcomp::cld() tratamento Ano.trend SE df lower.CL upper.CL .group
100 0.411 0.184 36 0.0385 0.784 1
N25 0.444 0.184 36 0.0710 0.817 1
F50 0.562 0.184 36 0.1894 0.935 1
N50 0.638 0.184 36 0.2656 1.011 1
F25 0.780 0.184 36 0.4076 1.153 1
Results are averaged over the levels of: Ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. (iv) AGB only outside (unplanted areas)
Preparing data:
Comparison only outside plots over the years
show code
AGB_outside0 <- trees_clean |>
dplyr::group_by(ano, tratamento, bloco, parcela, plantio)|>
dplyr::summarise(
AGB_sum = sum(AGB_ferez, na.rm=T),
.groups = "drop") |>
dplyr::mutate(
tratamento = ordered(tratamento, levels = ordem),
Parcela = paste(tratamento, bloco, parcela, sep="_"),
Ano = ano-2021) |>
dplyr::filter(plantio == "SP")View
show code
AGB_outside0 |>
ggplot(aes(x = factor(ano), y = log(AGB_sum), group = Parcela, colour = parcela)) +
geom_point() +
geom_line()+facet_grid(bloco~tratamento)+
theme_minimal()
it is a mess
Tratar os valores como (1|bloco)?
show code
AGB_outside <- trees_clean |>
dplyr::filter(plantio == "SP") |>
#dplyr::filter(ano == 2024) |>
dplyr::group_by(ano, tratamento, bloco)|>
dplyr::summarise(AGB_sum = sum(AGB_ferez, na.rm=T)/1000, .groups = "drop") |>
dplyr::mutate(
Tratamento = ordered(tratamento, levels = ordem),
Ano = ano-2021,
ANO = ordered(ano),
scaled_AGB = scale(AGB_sum)) Comparing models
show code
AGB_outside_m0 <- lme4::lmer(scaled_AGB ~ 0 + (1|bloco), data = AGB_outside) # best modelboundary (singular) fit: see help('isSingular')show code
AGB_outside_m1 <- lme4::lmer(scaled_AGB ~ ANO + (1|bloco), data = AGB_outside)
AGB_outside_m2 <- lme4::lmer(scaled_AGB ~ ANO + bloco + (1|bloco), data = AGB_outside) # best model
AGB_outside_m3 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + bloco + (1|bloco), data = AGB_outside)
AGB_outside_m4 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + bloco + (1|bloco), data = AGB_outside) Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
unable to evaluate scaled gradientWarning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Hessian is numerically singular: parameters are not uniquely determinedshow code
AGB_outside_m5 <- lme4::lmer(scaled_AGB ~ tratamento + ANO + (1|bloco), data = AGB_outside)
AGB_outside_m6 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1|bloco), data = AGB_outside)
AGB_outside_m7 <- lme4::lmer(scaled_AGB ~ tratamento*ANO + (1+Ano|bloco), data = AGB_outside) # isSingularboundary (singular) fit: see help('isSingular')show code
AIC(AGB_outside_m0, AGB_outside_m1, AGB_outside_m2, AGB_outside_m3,
AGB_outside_m4,AGB_outside_m5, AGB_outside_m6, AGB_outside_m7) df AIC
AGB_outside_m0 2 116.5024
AGB_outside_m1 4 117.0739
AGB_outside_m2 8 120.3459
AGB_outside_m3 14 128.0330
AGB_outside_m4 11 126.0844
AGB_outside_m5 7 122.8125
AGB_outside_m6 10 124.7610
AGB_outside_m7 12 126.6309Exploring residuals
show code
par(mfrow=c(1,2))
DHARMa::plotResiduals(AGB_outside_m1)
DHARMa::plotQQunif(AGB_outside_m0)
Marginal \(R^2\) / Conditional \(R^2\)
show code
MuMIn::r.squaredGLMM(AGB_outside_m1) R2m R2c
[1,] 0.1505529 0.1714753show code
0.1505529 /0.1714753[1] 0.87798590.1505529 /0.1714753 variance explained by fixed effects and by the entire model (around 0.8779859 % is explained by fixed effects)
show code
car::Anova(AGB_outside_m5)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: scaled_AGB
Chisq Df Pr(>Chisq)
tratamento 0.6275 3 0.89012
ANO 6.5922 1 0.01024 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1interaction is no significant at all (\(p-value \sim 0.85\))
show code
sjPlot::tab_model(AGB_outside_m5)| scaled AGB | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 0.20 | -0.42 – 0.83 | 0.514 |
| tratamento [F50] | -0.30 | -1.17 – 0.57 | 0.493 |
| tratamento [N25] | -0.23 | -1.10 – 0.64 | 0.597 |
| tratamento [N50] | -0.29 | -1.15 – 0.58 | 0.509 |
| ANO [linear] | 0.55 | 0.11 – 0.98 | 0.015 |
| Random Effects | |||
| σ2 | 0.91 | ||
| τ00 bloco | 0.01 | ||
| ICC | 0.01 | ||
| N bloco | 5 | ||
| Observations | 40 | ||
| Marginal R2 / Conditional R2 | 0.154 / 0.167 | ||
show code
broom.mixed::tidy(AGB_outside_m5) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 7 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 0.202 0.307 0.66
2 fixed <NA> tratamentoF50 -0.296 0.427 -0.694
3 fixed <NA> tratamentoN25 -0.228 0.427 -0.534
4 fixed <NA> tratamentoN50 -0.285 0.427 -0.667
5 fixed <NA> ANO.L 0.549 0.214 2.57
6 ran_pars bloco sd__(Intercept) 0.116 NA NA
7 ran_pars Residual sd__Observation 0.956 NA NA Assessing comparisons:
show code
emmeans::emmeans(AGB_outside_m6, pairwise ~ tratamento |ANO) |> multcomp::cld()ANO = 2022:
tratamento emmean SE df lower.CL upper.CL .group
F50 -0.6400 0.436 32 -1.5273 0.247 1
N25 -0.4221 0.436 32 -1.3094 0.465 1
F25 -0.3904 0.436 32 -1.2776 0.497 1
N50 -0.0995 0.436 32 -0.9868 0.788 1
ANO = 2024:
tratamento emmean SE df lower.CL upper.CL .group
N50 -0.0659 0.436 32 -0.9531 0.821 1
N25 0.3708 0.436 32 -0.5165 1.258 1
F50 0.4519 0.436 32 -0.4354 1.339 1
F25 0.7952 0.436 32 -0.0921 1.682 1
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 4 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Table S2
Only inside
show code
AGB_inside |>
dplyr::group_by(ano, tratamento) |>
dplyr::summarise(
Mean = round(mean(AGB_sum),2),
SE = round(sd(AGB_sum)/sqrt(dplyr::n()), 2),
.groups="drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)# A tibble: 2 × 6
ano `100` F25 F50 N25 N50
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 1.23 pm 0.15 1.32 pm 0.2 1.12 pm 0.26 0.85 pm 0.2 0.96 pm 0.12
2 2024 1.85 pm 0.18 2.5 pm 0.37 1.96 pm 0.46 1.52 pm 0.37 1.92 pm 0.24Increment per year:
show code
AGB_inside |>
dplyr::group_by(tratamento) |>
dplyr::summarise(
variacao_media = mean((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
desvio_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
erro_padrao = sd((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2) / sqrt(dplyr::n()),
min_variacao = min((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
max_variacao = max((AGB_sum[ano == 2024] - AGB_sum[ano == 2022])/2),
.groups="drop")|>
dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
dplyr::select(tratamento, slope)# A tibble: 5 × 2
tratamento slope
<chr> <chr>
1 100 0.31 pm 0.03
2 F25 0.59 pm 0.08
3 F50 0.42 pm 0.08
4 N25 0.33 pm 0.06
5 N50 0.48 pm 0.1 Total Aboveground biomass
show code
tabS1_agb_total <- AGB_data |>
dplyr::group_by(ano, tratamento)|>
dplyr::summarise(Mean = round(mean(AGB_Mgha),2),
SE = round(sd(AGB_Mgha)/sqrt(dplyr::n()), 2),
.groups="drop")
tabS1_agb_total |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = tratamento, values_from = values)# A tibble: 2 × 6
ano `100` F25 F50 N25 N50
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 21.35 pm 2.64 6.47 pm 0.92 10.03 pm 2.29 4.37 pm 1.07 9.06 pm 0.95
2 2024 32.06 pm 3.11 12.81 pm 2.12 18.12 pm 4.4 8.11 pm 1.64 17.39 pm 1.95Increment per year:
show code
AGB_data |>
dplyr::group_by(tratamento) %>%
dplyr::summarise(
variacao_media = mean((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
desvio_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
erro_padrao = sd((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2) / sqrt(dplyr::n()),
min_variacao = min((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
max_variacao = max((AGB_Mgha[ano == 2024] - AGB_Mgha[ano == 2022])/2),
.groups="drop")|>
dplyr::mutate(slope = paste(round(variacao_media,2), "pm", round(erro_padrao,2)))|>
dplyr::select(tratamento, slope)# A tibble: 5 × 2
tratamento slope
<chr> <chr>
1 100 5.35 pm 0.59
2 F25 3.17 pm 0.51
3 F50 4.05 pm 0.82
4 N25 1.87 pm 0.35
5 N50 4.17 pm 0.9 Figure Biomass
Percentuals
show code
agb_plot_data <-AGB_data0 |>
dplyr::group_by(ano, Tratamento, plantio, bloco)|>
dplyr::summarise(
AGB_Mgha = sum(AGB_Mgha),
.groups = "drop") |>
dplyr::group_by(ano, Tratamento, plantio)|>
dplyr::summarise(meanAGB = mean(AGB_Mgha),
sd = sd(AGB_Mgha),
n = dplyr::n(),
se = sd/sqrt(n),
.groups = "drop") |>
dplyr::mutate(
Ano = as.factor(ifelse(ano == 2022, "5", "7")),
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
))
max_year <- agb_plot_data |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(meanAGB), .groups="drop")
agb_plot_data |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(meanAGB = max(meanAGB),.groups="drop")|>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (meanAGB/max_year)*100, 0))# A tibble: 2 × 5
ano Tratamento meanAGB max_year porcentos
<dbl> <ord> <dbl> <dbl> <dbl>
1 2022 <NA> 48.9 48.9 100
2 2024 <NA> 83.2 83.2 100Plotting
show code
error_bars <- agb_plot_data |>
arrange(Ano, Tratamento, desc(plantio)) |>
group_by(Ano, Tratamento) |>
mutate(meanAGB_new = cumsum(meanAGB)) |>
ungroup()
porcentos <- data.frame(
Ano = rep(c("5", "7"), each = 5),
Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
meanAGB = c(8, 9, 12, 15, 27, 12, 17, 23, 25, 38),
label = c(c("17%","27%", "39%", "46%", "100%",
"21%", "34%", "52%", "53%", "100%"))
)
text_data <- data.frame(
Ano = rep(c("5", "7"), each = 5),
Tratamento = rep(unique(agb_plot_data$Tratamento), times = length(unique(agb_plot_data$Ano))),
meanAGB = c(10, 11, 14, 17, 29, 14, 19, 25, 27, 40),
label = c(c("a","a", "a", "a", "b",
"A", "AB", "B", "B", "C"))
)
fig_biomass<-agb_plot_data |>
ggplot(aes(x = Ano, y = meanAGB, fill = plantio)) +
geom_bar(stat = "identity", color = "black") +
geom_errorbar(data=error_bars,
aes(ymin = meanAGB_new - se, ymax = meanAGB_new + se),
width = .3) +
scale_fill_manual(
labels = c('Planted', 'Unplanted'),
values = c("#92D050", "white")
) +
labs(fill = "",
y = expression("Total aboveground biomass (Mg.ha"^{-1} ~")"),
x = "") +
facet_grid(.~Tratamento, switch = "both")+
geom_text(data = text_data, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
geom_text(data = porcentos, aes(x = Ano, y = meanAGB, label = label), size = 3, color = "black", inherit.aes = FALSE) +
scale_y_continuous(expand = c(0, 0), limits = c(0,43))+
ggthemes::theme_clean()+
theme(
#plot.margin = margin(b = 0, l = 0),
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_blank(),
axis.text.x = element_text(color="black", hjust = 1, size=10),
axis.text.y = element_text(color="black", size=9),
axis.title.y = element_text(color="black", size=11),
strip.text = element_text(color="black", hjust = 0.5, size=10),
panel.spacing = unit(0.1, "mm"),
strip.placement = "outside",
legend.text = element_text(size = 11),
legend.title = element_blank(),
legend.key.size = unit(3, "mm"),
legend.box.background = element_blank(),
legend.box.margin = margin(t = 5, r = 0, b = 0, l = 0, "mm"),
legend.background = element_blank(),
legend.margin = margin(0,0,0,0),
legend.position = c(0.18, 0.90)
)
fig_biomassWarning: Removed 2 rows containing missing values or values outside the scale range
(`geom_bar()`).
show code
# ggsave("fig_biomass.png", plot=fig_biomass, width=90, height = 90 ,unit="mm")II. Results - LiDAR
Canopy Height Metric (CHM_mean, CHM_SD, and CHM_CV), Canopy cover (Canopy_cover), Rugosity (Rugosity_SD, and Rugosity_CV), Leaf Area Index (LAI_mean, LAI_SD, and LAI_CV), and LAIunder (LAIunder_mean, LAIunder_SD, LAIunder_CV).
a) Canopy Height Model
2022
show code
CHM2022 <- lidar_wide |>
dplyr::filter(ano==2022)|>
with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2022)
DHARMa::plotQQunif(CHM2022)
2024
show code
CHM2024 <- lidar_wide |>
dplyr::filter(ano==2024)|>
with(lme4::lmer(CHM_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(CHM2024)
DHARMa::plotQQunif(CHM2024)
Together
show code
modeloCHM <- lidar_wide |>
with(lme4::lmer(CHM_mean ~ Tratamento*ano + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCHM)
DHARMa::plotQQunif(modeloCHM)
Comparing estimated means
show code
emmeans::emmeans(modeloCHM, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 1.89 0.437 10 0.919 2.87 1
Strips 25% 2.12 0.437 10 1.142 3.09 1
Strips 50% 2.93 0.437 10 1.956 3.90 1
Nuclei 50% 3.00 0.437 10 2.027 3.97 1
Planting 100% 4.91 0.437 10 3.941 5.89 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 2.76 0.437 10 1.788 3.73 1
Strips 25% 3.16 0.437 10 2.187 4.13 1
Strips 50% 4.43 0.437 10 3.454 5.40 2
Nuclei 50% 4.70 0.437 10 3.724 5.67 2
Planting 100% 6.99 0.437 10 6.018 7.96 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Comparing estimated slopes
show code
emmeans::emtrends(modeloCHM, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% 0.434 0.201 36 0.026 0.843 1
Strips 25% 0.523 0.201 36 0.114 0.931 1
Strips 50% 0.749 0.201 36 0.341 1.157 1
Nuclei 50% 0.849 0.201 36 0.440 1.257 1
Planting 100% 1.039 0.201 36 0.630 1.447 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Comparing fixed and random factors variance
show code
MuMIn::r.squaredGLMM(modeloCHM) R2m R2c
[1,] 0.7002626 0.8728069show code
0.7002626 /0.8728069[1] 0.802311Results:
The trend of variation around the mean CHM was the same over time among treatments (n.s.). The significant differences (\(p<0.001\)) observed between the 100% planting treatment and the other treatments in 2022 persisted in 2024.
show code
sjPlot::tab_model(modeloCHM)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -876.29 | -1700.80 – -51.79 | 0.038 |
| Tratamento [Strips 25%] | -178.32 | -1344.35 – 987.70 | 0.759 |
| Tratamento [Nuclei 50%] | -836.96 | -2002.98 – 329.07 | 0.154 |
| Tratamento [Strips 50%] | -635.54 | -1801.57 – 530.48 | 0.277 |
| Tratamento [Planting 100%] |
-1218.82 | -2384.85 – -52.80 | 0.041 |
| ano | 0.43 | 0.03 – 0.84 | 0.037 |
| Tratamento [Strips 25%] × ano |
0.09 | -0.49 – 0.66 | 0.758 |
| Tratamento [Nuclei 50%] × ano |
0.41 | -0.16 – 0.99 | 0.154 |
| Tratamento [Strips 50%] × ano |
0.31 | -0.26 – 0.89 | 0.276 |
| Tratamento [Planting 100%] × ano |
0.60 | 0.03 – 1.18 | 0.040 |
| Random Effects | |||
| σ2 | 0.41 | ||
| τ00 Bloco | 0.55 | ||
| ICC | 0.58 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.700 / 0.873 | ||
show code
broom.mixed::tidy(modeloCHM) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -876. 407. -2.15
2 fixed <NA> TratamentoStrips 25% -178. 576. -0.31
3 fixed <NA> TratamentoNuclei 50% -837. 576. -1.45
4 fixed <NA> TratamentoStrips 50% -636. 576. -1.10
5 fixed <NA> TratamentoPlanting 100% -1219. 576. -2.12
6 fixed <NA> ano 0.434 0.201 2.16
7 fixed <NA> TratamentoStrips 25%:ano 0.088 0.285 0.31
8 fixed <NA> TratamentoNuclei 50%:ano 0.414 0.285 1.46
9 fixed <NA> TratamentoStrips 50%:ano 0.315 0.285 1.11
10 fixed <NA> TratamentoPlanting 100%:ano 0.604 0.285 2.12
11 ran_pars Bloco sd__(Intercept) 0.742 NA NA
12 ran_pars Residual sd__Observation 0.637 NA NA Table S2
Values per year
show code
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(CHM_mean), 2),
SE = round(sd(CHM_mean)/sqrt(dplyr::n()), 2)) |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)`summarise()` has grouped output by 'ano'. You can override using the `.groups`
argument.# A tibble: 2 × 6
# Groups: ano [2]
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 1.89 pm 0.3 2.12 pm 0.41 3 pm 0.26 2.93 pm 0.46 4.91 pm 0.33
2 2024 2.76 pm 0.48 3.16 pm 0.56 4.7 pm 0.28 4.43 pm 0.69 6.99 pm 0.4 “Slope” (increment/year)
show code
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(CHM_mean, na.rm = TRUE),
se = sd(CHM_mean, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.43 pm 0.28 0.52 pm 0.35 0.85 pm 0.19 0.75 pm 0.42 1.04 pm 0.26 b) Canopy cover
show code
lidar_wide |>
dplyr::group_by(ano)|>
dplyr::summarise(mean(Canopy_cover),
min(Canopy_cover),
max(Canopy_cover),
sd(Canopy_cover))# A tibble: 2 × 5
ano `mean(Canopy_cover)` `min(Canopy_cover)` `max(Canopy_cover)`
<dbl> <dbl> <dbl> <dbl>
1 2022 38.4 5.29 99.8
2 2024 67.6 10.7 99.8
# ℹ 1 more variable: `sd(Canopy_cover)` <dbl>2022
show code
CC2022 <- lidar_wide |>
dplyr::filter(ano==2022)|>
with( glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1|Bloco)) )Warning in glmmTMB::glmmTMB(Canopy_cover ~ Tratamento + (1 | Bloco)): use of
the 'data' argument is recommendedshow code
CC2022a <- lidar_wide |>
dplyr::filter(ano==2022)|>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommendedshow code
AIC(CC2022, CC2022a) df AIC
CC2022 7 228.16951
CC2022a 7 -5.60719show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2022a)
DHARMa::plotQQunif(CC2022a)
DHARMa::plotResiduals(CC2022)
DHARMa::plotQQunif(CC2022)
2024
show code
CC2024 <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(lme4::lmer(Canopy_cover ~ Tratamento + (1|Bloco) ) )
CC2024a <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento + (1 | Bloco), : use
of the 'data' argument is recommendedshow code
AIC(CC2024, CC2024a) df AIC
CC2024 7 212.04388
CC2024a 7 -16.11127show code
par(mfrow=c(2,2))
DHARMa::plotResiduals(CC2024a)
DHARMa::plotQQunif(CC2024a)
DHARMa::plotResiduals(CC2024)
DHARMa::plotQQunif(CC2024)
models
show code
modeloCC <- lidar_wide |>
with(lme4::lmer(Canopy_cover ~ Tratamento*ano + (1|Bloco) ) )
modeloCCa <- lidar_wide |>
with(glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento*ano + (1 | Bloco),
family = glmmTMB::beta_family()))Warning in glmmTMB::glmmTMB(Canopy_cover/100 ~ Tratamento * ano + (1 | Bloco),
: use of the 'data' argument is recommendedWarning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluation
Warning in (function (start, objective, gradient = NULL, hessian = NULL, :
NA/NaN function evaluationshow code
par(mfrow=c(2,2))
DHARMa::plotResiduals(modeloCC)
DHARMa::plotQQunif(modeloCC)
DHARMa::plotResiduals(modeloCCa)
DHARMa::plotQQunif(modeloCCa)
1 outlier
Estimated means
show code
emmeans::emmeans(modeloCC, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Strips 25% 23.7 11.5 20.8 -0.25 47.6 1
Nuclei 50% 26.5 11.5 20.8 2.56 50.4 1
Strips 50% 31.3 11.5 20.8 7.42 55.3 1
Nuclei 25% 31.6 11.5 20.8 7.70 55.5 1
Planting 100% 79.1 11.5 20.8 55.14 103.0 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Strips 25% 44.7 11.5 20.8 20.78 68.6 1
Nuclei 25% 56.5 11.5 20.8 32.62 80.5 1
Strips 50% 69.8 11.5 20.8 45.88 93.7 12
Nuclei 50% 70.9 11.5 20.8 47.02 94.9 12
Planting 100% 96.2 11.5 20.8 72.29 120.1 2
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Estimated slopes
show code
emmeans::emtrends(modeloCC, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Planting 100% 8.57 6.71 36 -5.02 22.2 1
Strips 25% 10.52 6.71 36 -3.08 24.1 1
Nuclei 25% 12.46 6.71 36 -1.14 26.1 1
Strips 50% 19.23 6.71 36 5.63 32.8 1
Nuclei 50% 22.23 6.71 36 8.63 35.8 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Comparing random and fixed effects
show code
MuMIn::r.squaredGLMM(modeloCC) R2m R2c
[1,] 0.469023 0.6387404show code
0.469023 / 0.6387404[1] 0.7342936show code
sjPlot::tab_model(modeloCC)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -25156.96 | -52617.32 – 2303.39 | 0.071 |
| Tratamento [Strips 25%] | 3917.03 | -34917.77 – 42751.83 | 0.839 |
| Tratamento [Nuclei 50%] | -19762.19 | -58596.99 – 19072.61 | 0.309 |
| Tratamento [Strips 50%] | -13690.08 | -52524.88 – 25144.72 | 0.480 |
| Tratamento [Planting 100%] |
7898.75 | -30936.05 – 46733.55 | 0.683 |
| ano | 12.46 | -1.12 – 26.03 | 0.071 |
| Tratamento [Strips 25%] × ano |
-1.94 | -21.14 – 17.26 | 0.839 |
| Tratamento [Nuclei 50%] × ano |
9.77 | -9.43 – 28.97 | 0.309 |
| Tratamento [Strips 50%] × ano |
6.77 | -12.43 – 25.97 | 0.480 |
| Tratamento [Planting 100%] × ano |
-3.88 | -23.08 – 15.31 | 0.684 |
| Random Effects | |||
| σ2 | 449.60 | ||
| τ00 Bloco | 211.22 | ||
| ICC | 0.32 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.469 / 0.639 | ||
show code
broom.mixed::tidy(modeloCC) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -25157. 13565. -1.86
2 fixed <NA> TratamentoStrips 25% 3917. 19183. 0.204
3 fixed <NA> TratamentoNuclei 50% -19762. 19183. -1.03
4 fixed <NA> TratamentoStrips 50% -13690. 19183. -0.714
5 fixed <NA> TratamentoPlanting 100% 7899. 19183. 0.412
6 fixed <NA> ano 12.5 6.70 1.86
7 fixed <NA> TratamentoStrips 25%:ano -1.94 9.48 -0.205
8 fixed <NA> TratamentoNuclei 50%:ano 9.77 9.48 1.03
9 fixed <NA> TratamentoStrips 50%:ano 6.77 9.48 0.714
10 fixed <NA> TratamentoPlanting 100%:ano -3.88 9.48 -0.409
11 ran_pars Bloco sd__(Intercept) 14.5 NA NA
12 ran_pars Residual sd__Observation 21.2 NA NA Table S2
per year
show code
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(Canopy_cover), 2),
SE = round((sd(Canopy_cover)/sqrt(dplyr::n())), 2),
.groups = "drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 31.62 pm 17.23 23.67 pm 4.99 26.48 pm 1.81 31.34 pm 7.74 79.06 pm 10.54
2 2024 56.54 pm 18.23 44.7 pm 13.14 70.94 pm 11.3 69.8 pm 13.76 96.2 pm 1.89 “Slope” (increment/year)
show code
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(Canopy_cover, na.rm = TRUE),
se = sd(Canopy_cover, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 12.46 pm 12.54 10.52 pm 7.03 22.23 pm 5.72 19.23 pm 7.89 8.57 pm 5.35 c) variance of Canopy cover
Performing a pairwise comparison of variances is not as straightforward as comparing means, but you can use specific statistical techniques to compare variances between treatment groups. Here’s how you can do this in R:
show code
pairwise_levene <- function(data, response, group) {
treatment_pairs <- combn(unique(data[[group]]), 2, simplify = FALSE)
p_values <- vector("list", length(treatment_pairs))
for (i in seq_along(treatment_pairs)) { #i =1
subset_data <- data[data[[group]] %in% treatment_pairs[[i]], ]
test_result <- leveneTest(subset_data[[response]] ~ subset_data[[group]])
p_values[[i]] <- test_result[1, "Pr(>F)"]
}
comparison <- sapply(treatment_pairs, function(x) paste(x, collapse = " vs "))
results <- data.frame(Comparison = comparison, P_Value = unlist(p_values))
results$Adjusted_P_Value <- p.adjust(results$P_Value, method = "bonferroni")
return(results)
}
pairwise_levene(lidar2022_wide , "Canopy_cover", "Tratamento")Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor. Comparison P_Value Adjusted_P_Value
1 25_faixa vs 25_tabuleiro 0.41911844 1.0000000
2 25_faixa vs 50_faixa 0.60465350 1.0000000
3 25_faixa vs 50_tabuleiro 0.52704716 1.0000000
4 25_faixa vs 100_controle 0.21179066 1.0000000
5 25_tabuleiro vs 50_faixa 0.56742080 1.0000000
6 25_tabuleiro vs 50_tabuleiro 0.31075886 1.0000000
7 25_tabuleiro vs 100_controle 0.89320378 1.0000000
8 50_faixa vs 50_tabuleiro 0.26004211 1.0000000
9 50_faixa vs 100_controle 0.43307531 1.0000000
10 50_tabuleiro vs 100_controle 0.07724271 0.7724271show code
pairwise_levene(lidar2024_wide , "Canopy_cover", "Tratamento")Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor.
Warning in leveneTest.default(y = y, group = group, ...): group coerced to
factor. Comparison P_Value Adjusted_P_Value
1 25_faixa vs 25_tabuleiro 0.47523130 1.0000000
2 25_faixa vs 50_faixa 0.73908327 1.0000000
3 25_faixa vs 50_tabuleiro 0.87724043 1.0000000
4 25_faixa vs 100_controle 0.28751774 1.0000000
5 25_tabuleiro vs 50_faixa 0.65306786 1.0000000
6 25_tabuleiro vs 50_tabuleiro 0.49847386 1.0000000
7 25_tabuleiro vs 100_controle 0.08104640 0.8104640
8 50_faixa vs 50_tabuleiro 0.81722803 1.0000000
9 50_faixa vs 100_controle 0.11052807 1.0000000
10 50_tabuleiro vs 100_controle 0.09110862 0.9110862Table S2 (Rugosity)
show code
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(Rugosity_SD),2),
SE = round((sd(Rugosity_SD)/sqrt(dplyr::n())), 2),
.groups="drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 2.08 pm 0.22 2.38 pm 0.17 2.73 pm 0.31 2.64 pm 0.19 2.57 pm 0.06
2 2024 2.77 pm 0.32 3.15 pm 0.18 3.41 pm 0.36 3.3 pm 0.17 2.79 pm 0.09 Increment/year
show code
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(Rugosity_SD, na.rm = TRUE),
se = sd(Rugosity_SD, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.34 pm 0.19 0.39 pm 0.12 0.34 pm 0.24 0.33 pm 0.13 0.11 pm 0.05 d) LAI Mean
2022
show code
LAI2022 <- lidar_wide |>
dplyr::filter(ano==2022) |>
with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2022)
DHARMa::plotQQunif(LAI2022)
2024
show code
LAI2024 <- lidar_wide |>
dplyr::filter(ano==2024) |>
with(lme4::lmer(LAI_mean ~ Tratamento + (1|Bloco) ) )
par(mfrow=c(1,2))
DHARMa::plotResiduals(LAI2024)
DHARMa::plotQQunif(LAI2024)
Full model
show code
modeloLAI <- lidar_wide |>
with(lme4::lmer(LAI_mean ~ Tratamento * ano + (1 | Bloco)))
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloLAI)
DHARMa::plotQQunif(modeloLAI)
Pairwise comparison
show code
emmeans::emmeans(modeloLAI, pairwise ~ Tratamento | ano) |> multcomp::cld()ano = 2022:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 0.447 0.159 13.2 0.105 0.789 1
Strips 25% 0.591 0.159 13.2 0.249 0.933 1
Nuclei 50% 0.781 0.159 13.2 0.439 1.123 1
Strips 50% 0.828 0.159 13.2 0.485 1.170 1
Planting 100% 1.565 0.159 13.2 1.223 1.907 2
ano = 2024:
Tratamento emmean SE df lower.CL upper.CL .group
Nuclei 25% 0.900 0.159 13.2 0.558 1.242 1
Strips 25% 1.047 0.159 13.2 0.705 1.389 12
Strips 50% 1.353 0.159 13.2 1.011 1.695 12
Nuclei 50% 1.383 0.159 13.2 1.041 1.725 2
Planting 100% 2.408 0.159 13.2 2.066 2.750 3
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Trend:
show code
emmeans::emtrends(modeloLAI, pairwise ~ Tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() Tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% 0.226 0.0812 36 0.0617 0.391 1
Strips 25% 0.228 0.0812 36 0.0633 0.393 1
Strips 50% 0.263 0.0812 36 0.0980 0.428 1
Nuclei 50% 0.301 0.0812 36 0.1360 0.466 1
Planting 100% 0.421 0.0812 36 0.2565 0.586 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. Fixed vs random
show code
MuMIn::r.squaredGLMM(modeloLAI) R2m R2c
[1,] 0.7064328 0.8459851show code
0.7064328 /0.8459851[1] 0.8350417show code
sjPlot::tab_model(modeloLAI)| Dependent variable | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | -457.36 | -790.02 – -124.71 | 0.008 |
| Tratamento [Strips 25%] | -3.08 | -473.52 – 467.36 | 0.989 |
| Tratamento [Nuclei 50%] | -150.02 | -620.46 – 320.42 | 0.522 |
| Tratamento [Strips 50%] | -73.12 | -543.56 – 397.32 | 0.755 |
| Tratamento [Planting 100%] |
-392.88 | -863.32 – 77.56 | 0.099 |
| ano | 0.23 | 0.06 – 0.39 | 0.008 |
| Tratamento [Strips 25%] × ano |
0.00 | -0.23 – 0.23 | 0.989 |
| Tratamento [Nuclei 50%] × ano |
0.07 | -0.16 – 0.31 | 0.521 |
| Tratamento [Strips 50%] × ano |
0.04 | -0.20 – 0.27 | 0.753 |
| Tratamento [Planting 100%] × ano |
0.19 | -0.04 – 0.43 | 0.098 |
| Random Effects | |||
| σ2 | 0.07 | ||
| τ00 Bloco | 0.06 | ||
| ICC | 0.48 | ||
| N Bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.706 / 0.846 | ||
show code
broom.mixed::tidy(modeloLAI) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), ~ round(.x, 3)))# A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) -457. 164. -2.78
2 fixed <NA> TratamentoStrips 25% -3.08 232. -0.013
3 fixed <NA> TratamentoNuclei 50% -150. 232. -0.646
4 fixed <NA> TratamentoStrips 50% -73.1 232. -0.315
5 fixed <NA> TratamentoPlanting 100% -393. 232. -1.69
6 fixed <NA> ano 0.226 0.081 2.79
7 fixed <NA> TratamentoStrips 25%:ano 0.002 0.115 0.014
8 fixed <NA> TratamentoNuclei 50%:ano 0.074 0.115 0.647
9 fixed <NA> TratamentoStrips 50%:ano 0.036 0.115 0.316
10 fixed <NA> TratamentoPlanting 100%:ano 0.195 0.115 1.70
11 ran_pars Bloco sd__(Intercept) 0.245 NA NA
12 ran_pars Residual sd__Observation 0.257 NA NA Table S2
show code
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = round(mean(LAI_mean),2),
SE = round((sd(LAI_mean)/sqrt(dplyr::n())), 2),
.groups = "drop") |>
dplyr::mutate(values = paste(Mean, "pm", SE)) |>
dplyr::select(-c(Mean, SE)) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 2 × 6
ano `Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<dbl> <chr> <chr> <chr> <chr> <chr>
1 2022 0.45 pm 0.09 0.59 pm 0.12 0.78 pm 0.08 0.83 pm 0.16 1.57 pm 0.17
2 2024 0.9 pm 0.15 1.05 pm 0.2 1.38 pm 0.16 1.35 pm 0.19 2.41 pm 0.21 Increment/year (slope)
show code
lidar_wide |>
dplyr::group_by(Tratamento, ano) |>
dplyr::summarise(
media = mean(LAI_mean, na.rm = TRUE),
se = sd(LAI_mean, na.rm = TRUE) / sqrt(dplyr::n()),
.groups = "drop"
) |>
tidyr::pivot_wider(
names_from = ano,
values_from = c(media, se),
names_sep = "_"
) |>
dplyr::mutate(
variacao_media = round((media_2024 - media_2022) / 2, 2),
SE_variacao = round(sqrt(se_2022^2 + se_2024^2) / 2, 2)) |>
dplyr::mutate(values = paste(variacao_media, "pm", SE_variacao)) |>
dplyr::select(Tratamento, values) |>
tidyr::pivot_wider(names_from = Tratamento, values_from = values)# A tibble: 1 × 5
`Nuclei 25%` `Strips 25%` `Nuclei 50%` `Strips 50%` `Planting 100%`
<chr> <chr> <chr> <chr> <chr>
1 0.23 pm 0.09 0.23 pm 0.12 0.3 pm 0.09 0.26 pm 0.12 0.42 pm 0.14 e) Figure - LiDAR
(i) Canopy height model (CHM)
show code
max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 1.89 4.91 39
2 2022 Strips 25% 2.12 4.91 43
3 2022 Nuclei 50% 3.00 4.91 61
4 2022 Strips 50% 2.93 4.91 60
5 2022 Planting 100% 4.91 4.91 100
6 2024 Nuclei 25% 2.76 6.99 39
7 2024 Strips 25% 3.16 6.99 45
8 2024 Nuclei 50% 4.70 6.99 67
9 2024 Strips 50% 4.43 6.99 63
10 2024 Planting 100% 6.99 6.99 100Plot
show code
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," a",
"A "," b",
"A "," b",
"B "," c"),
y = c(2.8, 3.9,
3.1, 4.5,
3.9, 5.7,
3.9, 5.7,
5.9, 7.8) )
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("39% "," 39%",
"43% "," 45%",
"61% "," 67%",
"60% "," 63%",
"100% "," 100%"),
y = c(2.2, 3.1,
2.5, 3.7,
3.3, 5,
3.3, 5,
5.3, 7.1)
)
p1 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(CHM_mean),
SD = sd(CHM_mean),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 9), breaks = seq(0, 8, by = 1))+
labs(fill = "",
y = "Canopy height model (m)",
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=9),
axis.title.y = element_text(color="black", size=10),
legend.position = c(0.15, 0.9),
legend.text = element_text(size = 10),
legend.key.size = unit(3, "mm"),
plot.margin = margin(t=5, 0, 0, 0, "mm"))
p1
(ii) Canopy cover
show code
max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 31.6 79.1 40
2 2022 Strips 25% 23.7 79.1 30
3 2022 Nuclei 50% 26.5 79.1 33
4 2022 Strips 50% 31.3 79.1 40
5 2022 Planting 100% 79.1 79.1 100
6 2024 Nuclei 25% 56.5 96.2 59
7 2024 Strips 25% 44.7 96.2 46
8 2024 Nuclei 50% 70.9 96.2 74
9 2024 Strips 50% 69.8 96.2 73
10 2024 Planting 100% 96.2 96.2 100plot
show code
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," a",
"A "," ab",
"A "," ab",
"B "," b"),
y = c(52, 77,
40, 65,
40, 85,
45, 90,
97, 109)
)
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("40% "," 59%",
"30% "," 46%",
"33% ", " 74%",
"40% "," 73%",
"100% "," 100%"),
y = c(42, 67,
30, 55,
30, 75,
35, 80,
87, 100)
)
p2 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(Canopy_cover),
SD = sd(Canopy_cover),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits=c(0, 120), breaks = seq(0, 100, by = 20))+
labs(fill = "",
y = "Canopy cover (%)",
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = 'none')
p2
(iii) LAI
show code
max_year <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
dplyr::group_by(ano)|>
dplyr::summarise(max_year = max(Mean),.groups="drop")
lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),.groups="drop") |>
dplyr::right_join(max_year, by="ano") |>
dplyr::mutate(porcentos = round( (Mean/max_year)*100, 0))# A tibble: 10 × 5
ano Tratamento Mean max_year porcentos
<dbl> <fct> <dbl> <dbl> <dbl>
1 2022 Nuclei 25% 0.447 1.57 29
2 2022 Strips 25% 0.591 1.57 38
3 2022 Nuclei 50% 0.781 1.57 50
4 2022 Strips 50% 0.828 1.57 53
5 2022 Planting 100% 1.57 1.57 100
6 2024 Nuclei 25% 0.900 2.41 37
7 2024 Strips 25% 1.05 2.41 43
8 2024 Nuclei 50% 1.38 2.41 57
9 2024 Strips 50% 1.35 2.41 56
10 2024 Planting 100% 2.41 2.41 100Plot
show code
letras <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("A "," a",
"A "," ab",
"A ", " b",
"A "," ab",
"B "," c"),
y = c( 0.85, 1.4,
1.05, 1.55,
1.15, 1.82,
1.25, 1.82,
2.05, 2.9)
) |>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
porcentos <- data.frame(
Tratamento = rep(levels(as.factor(lidar_wide$Tratamento)),each=2),
ano = c(2022, 2024),
label = c("29% "," 37%",
"38% "," 43%",
"50% "," 57%",
"53% "," 56%",
"100% "," 100%"),
y = c( 0.55, 1.1,
0.75, 1.25,
0.85, 1.52,
0.95, 1.52,
1.75, 2.65)
)|>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
)
)
p3 <- lidar_wide |>
dplyr::group_by(ano, Tratamento)|>
dplyr::summarise(Mean = mean(LAI_mean),
SD = sd(LAI_mean),
SE = SD/sqrt(dplyr::n()),
.groups="drop") |>
dplyr::mutate(
Tratamento = forcats::fct_recode(Tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"))|>
ggplot(aes(x = Tratamento, y = Mean, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = Mean - SE, ymax = Mean + SE), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = Tratamento, y = y, label = label),size = 3.5, color="black", vjust = -0.5) +
geom_text(data = porcentos, aes(x = Tratamento, y = y, label = label),size = 3, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0,3.2), breaks = seq(0, 3, by=1))+
labs(fill = "",
y = expression("LAI (m"^{2} ~ "/ m"^{2}~ ")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_text(color="black", hjust = 0.5, size=10),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = 'none')
p3
(iv) All together
show code
# 90 mm x 180 mm
fig_lidar <- cowplot::plot_grid(p1,p2,p3, ncol = 1,align = "v", axis = "l",
labels = c("(a)", "(b)", "(c)"), label_size = 10, label_x = 0.01, label_y = 1)
fig_lidar
show code
# ggsave("fig_lidar.png", plot=fig_lidar, width=90, height = 210 ,unit="mm")capital letters: 5 years
small letters: 7 years
III. Results - Costs
A) Cost per unit (Biomass and Canopy coverage)
Table of Costs (in R$)
show code
costs <- data.frame(
tratamento = c("Planting 100%", "Strips 50%", "Nuclei 50%", "Strips 25%", "Nuclei 25%"),
labor_machinery = c(7860.00, 6407.90, 7563.90, 5131.05, 6395.23),
inputs = c(7032.50, 4096.25, 4096.25, 2628.13, 2628.13)) |>
dplyr::mutate(total = labor_machinery + inputs)Using raw values
show code
cover_raw <- lidar_wide |>
dplyr::select(ano, Bloco, Tratamento, Canopy_cover) |>
dplyr::rename(bloco = Bloco, tratamento = Tratamento) |>
dplyr::mutate(tratamento = ordered(tratamento,
levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
bloco = ordered(bloco, levels = c("BI","BII", "BIII", "BIV", "BV")),
bloco = as.character(as.numeric(bloco)))|>
dplyr::left_join(costs) |>
dplyr::mutate(totalUS = (total/5)/5.8)Joining with `by = join_by(tratamento)`show code
cost_raw <- AGB_data |>
dplyr::select(ano, bloco, tratamento, AGB_Mgha) |>
dplyr::mutate(
tratamento = dplyr::case_when(
tratamento == "100" ~ "Planting 100%",
tratamento == "F25" ~ "Strips 25%",
tratamento == "F50" ~ "Strips 50%",
tratamento == "N25" ~ "Nuclei 25%",
tratamento == "N50" ~ "Nuclei 50%"
)
) |>
dplyr::right_join(cover_raw) |>
dplyr::mutate(agb_cost = totalUS/AGB_Mgha,
cc_cost = totalUS/Canopy_cover,
ano = ano-2017) |>
dplyr::select(ano:tratamento, agb_cost, cc_cost)Joining with `by = join_by(ano, bloco, tratamento)`(i) Tests
show code
modeloCC_cost <- lme4::lmer(log(cc_cost) ~ tratamento*ano + (1|bloco), data=cost_raw)
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloCC_cost)
DHARMa::plotQQunif(modeloCC_cost)
show code
emmeans::emmeans(modeloCC_cost, pairwise ~ tratamento | ano) |> multcomp::cld()ano = 5:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 1.91 0.256 18 1.37 2.45 1
Strips 25% 2.50 0.256 18 1.96 3.03 12
Strips 50% 2.56 0.256 18 2.02 3.10 12
Nuclei 50% 2.73 0.256 18 2.19 3.27 12
Nuclei 25% 2.77 0.256 18 2.24 3.31 2
ano = 7:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 1.68 0.256 18 1.14 2.21 1
Strips 50% 1.75 0.256 18 1.21 2.29 1
Nuclei 50% 1.79 0.256 18 1.25 2.32 1
Strips 25% 1.92 0.256 18 1.38 2.46 1
Nuclei 25% 1.99 0.256 18 1.45 2.53 1
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
sjPlot::tab_model(modeloCC_cost)| log(cc cost) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 4.73 | 2.93 – 6.54 | <0.001 |
| tratamento [Nuclei 50%] | 0.36 | -2.15 – 2.86 | 0.775 |
| tratamento [Planting 100%] |
-2.23 | -4.74 – 0.27 | 0.079 |
| tratamento [Strips 25%] | -0.81 | -3.31 – 1.70 | 0.519 |
| tratamento [Strips 50%] | -0.15 | -2.66 – 2.36 | 0.905 |
| ano | -0.39 | -0.68 – -0.10 | 0.010 |
| tratamento [Nuclei 50%] × ano |
-0.08 | -0.49 – 0.33 | 0.697 |
| tratamento [Planting 100%] × ano |
0.27 | -0.14 – 0.69 | 0.186 |
| tratamento [Strips 25%] × ano |
0.11 | -0.31 – 0.52 | 0.607 |
| tratamento [Strips 50%] × ano |
-0.01 | -0.42 – 0.40 | 0.952 |
| Random Effects | |||
| σ2 | 0.21 | ||
| τ00 bloco | 0.12 | ||
| ICC | 0.37 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.341 / 0.583 | ||
show code
broom.mixed::tidy(modeloCC_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) # A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 4.73 0.89 5.32
2 fixed <NA> tratamentoNuclei 50% 0.356 1.24 0.287
3 fixed <NA> tratamentoPlanting 100% -2.23 1.24 -1.80
4 fixed <NA> tratamentoStrips 25% -0.806 1.24 -0.651
5 fixed <NA> tratamentoStrips 50% -0.149 1.24 -0.12
6 fixed <NA> ano -0.392 0.144 -2.72
7 fixed <NA> tratamentoNuclei 50%:ano -0.08 0.204 -0.392
8 fixed <NA> tratamentoPlanting 100%:ano 0.274 0.204 1.35
9 fixed <NA> tratamentoStrips 25%:ano 0.106 0.204 0.518
10 fixed <NA> tratamentoStrips 50%:ano -0.012 0.204 -0.061
11 ran_pars bloco sd__(Intercept) 0.347 NA NA
12 ran_pars Residual sd__Observation 0.455 NA NA show code
modeloAGB_cost <- lme4::lmer(log(agb_cost) ~ tratamento*ano + (1|bloco), data=cost_raw )
par(mfrow=c(1,2))
DHARMa::plotResiduals(modeloAGB_cost)
DHARMa::plotQQunif(modeloAGB_cost)
show code
emmeans::emmeans(modeloAGB_cost, pairwise ~ tratamento | ano) |> multcomp::cld()ano = 5:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 3.21 0.217 17.5 2.75 3.67 1
Strips 50% 3.68 0.217 17.5 3.22 4.14 1
Strips 25% 3.76 0.217 17.5 3.30 4.22 1
Nuclei 50% 3.81 0.217 17.5 3.36 4.27 12
Nuclei 25% 4.50 0.217 17.5 4.04 4.95 2
ano = 7:
tratamento emmean SE df lower.CL upper.CL .group
Planting 100% 2.79 0.217 17.5 2.33 3.25 1
Strips 25% 3.08 0.217 17.5 2.63 3.54 1
Strips 50% 3.12 0.217 17.5 2.67 3.58 12
Nuclei 50% 3.17 0.217 17.5 2.71 3.62 12
Nuclei 25% 3.80 0.217 17.5 3.35 4.26 2
Degrees-of-freedom method: kenward-roger
Results are given on the log (not the response) scale.
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. And the slopes?
show code
emmeans::emtrends(modeloCC_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 50% -0.472 0.144 36 -0.764 -0.18000 1
Strips 50% -0.405 0.144 36 -0.697 -0.11250 1
Nuclei 25% -0.392 0.144 36 -0.684 -0.10005 1
Strips 25% -0.287 0.144 36 -0.579 0.00551 1
Planting 100% -0.118 0.144 36 -0.410 0.17419 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
emmeans::emtrends(modeloAGB_cost, pairwise ~ tratamento, var = "ano", poly.trend = TRUE)|> multcomp::cld() tratamento ano.trend SE df lower.CL upper.CL .group
Nuclei 25% -0.346 0.121 36 -0.592 -0.1006 1
Strips 25% -0.338 0.121 36 -0.584 -0.0928 1
Nuclei 50% -0.324 0.121 36 -0.570 -0.0782 1
Strips 50% -0.280 0.121 36 -0.525 -0.0340 1
Planting 100% -0.210 0.121 36 -0.456 0.0355 1
Results are averaged over the levels of: ano
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
P value adjustment: tukey method for comparing a family of 5 estimates
significance level used: alpha = 0.05
NOTE: If two or more means share the same grouping symbol,
then we cannot show them to be different.
But we also did not show them to be the same. show code
sjPlot::tab_model(modeloAGB_cost)| log(agb cost) | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 6.23 | 4.71 – 7.74 | <0.001 |
| tratamento [Nuclei 50%] | -0.79 | -2.90 – 1.32 | 0.452 |
| tratamento [Planting 100%] |
-1.96 | -4.07 – 0.15 | 0.067 |
| tratamento [Strips 25%] | -0.77 | -2.88 – 1.34 | 0.463 |
| tratamento [Strips 50%] | -1.15 | -3.26 – 0.96 | 0.278 |
| ano | -0.35 | -0.59 – -0.10 | 0.007 |
| tratamento [Nuclei 50%] × ano |
0.02 | -0.32 – 0.37 | 0.897 |
| tratamento [Planting 100%] × ano |
0.14 | -0.21 – 0.48 | 0.432 |
| tratamento [Strips 25%] × ano |
0.01 | -0.34 – 0.35 | 0.964 |
| tratamento [Strips 50%] × ano |
0.07 | -0.28 – 0.41 | 0.700 |
| Random Effects | |||
| σ2 | 0.15 | ||
| τ00 bloco | 0.09 | ||
| ICC | 0.38 | ||
| N bloco | 5 | ||
| Observations | 50 | ||
| Marginal R2 / Conditional R2 | 0.498 / 0.688 | ||
show code
broom.mixed::tidy(modeloAGB_cost) |> dplyr::mutate(dplyr::across(dplyr::where(is.numeric), round, 3)) # A tibble: 12 × 6
effect group term estimate std.error statistic
<chr> <chr> <chr> <dbl> <dbl> <dbl>
1 fixed <NA> (Intercept) 6.23 0.749 8.31
2 fixed <NA> tratamentoNuclei 50% -0.792 1.04 -0.76
3 fixed <NA> tratamentoPlanting 100% -1.96 1.04 -1.88
4 fixed <NA> tratamentoStrips 25% -0.773 1.04 -0.742
5 fixed <NA> tratamentoStrips 50% -1.15 1.04 -1.1
6 fixed <NA> ano -0.346 0.121 -2.86
7 fixed <NA> tratamentoNuclei 50%:ano 0.022 0.171 0.13
8 fixed <NA> tratamentoPlanting 100%:ano 0.136 0.171 0.794
9 fixed <NA> tratamentoStrips 25%:ano 0.008 0.171 0.045
10 fixed <NA> tratamentoStrips 50%:ano 0.067 0.171 0.389
11 ran_pars bloco sd__(Intercept) 0.298 NA NA
12 ran_pars Residual sd__Observation 0.383 NA NA (ii) Figure Cost-effectiveness
show code
fig_data_cost_raw <- cost_raw |>
dplyr::group_by(ano, tratamento)|>
dplyr::summarise(
cost_cc = mean(cc_cost),
SE_cc = sd(cc_cost)/sqrt(dplyr::n()),
cost_agb = mean(agb_cost),
SE_agb = sd(agb_cost)/sqrt(dplyr::n()),
.groups="drop")|>
dplyr::mutate(tratamento = ordered(tratamento, levels = c("Nuclei 25%", "Strips 25%", "Nuclei 50%", "Strips 50%", "Planting 100%")),
tratamento = forcats::fct_recode(tratamento,
"Nuclei\n25%" = "Nuclei 25%",
"Strips\n25%" = "Strips 25%",
"Nuclei\n50%" = "Nuclei 50%",
"Strips\n50%" = "Strips 50%",
"Planting\n100%" = "Planting 100%"
))show code
letras <- data.frame(
tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
ano = c(5, 7),
label = c("B "," a",
"AB "," a",
"AB "," a",
"AB "," a",
"A "," a"),
y = c(33, 16,
18, 10,
18, 8,
18, 10,
10, 6))
fig_cost1 <- fig_data_cost_raw |>
ggplot(aes(x = tratamento, y = cost_cc, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = cost_cc - SE_cc, ymax = cost_cc + SE_cc), width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
scale_y_continuous(expand=c(0,0), limits=c(0,38), breaks = c(0, 10,20, 30))+
labs(fill = "",
y = expression("CE Canopy cover (" * "\u0024" * ".%"^{-1} * ")"),
x = "") +
theme_classic()+
theme(
axis.text.x = element_blank(),
axis.text.y = element_text(color="black", size=10),
axis.title.y = element_text(color="black", size=10),
plot.margin = margin(t=5, 0, 0, 0, "mm"),
legend.position = c(0.85, 0.95))
fig_cost1
show code
letras <- data.frame(
tratamento = rep(levels(as.factor(fig_data_cost_raw$tratamento)), each=2),
ano = c(5, 7),
label = c("B "," b",
"A "," a",
"AB "," ab",
"A "," ab",
"A "," a"),
y = c(210, 100, 53, 25, 53, 30, 53, 33, 30, 20)
)
fig_cost2 <- fig_data_cost_raw |>
ggplot(aes(x = tratamento, y = cost_agb, fill = as.factor(ano))) +
geom_bar(position = position_dodge2(preserve = "single"),
stat = "identity", color = "black") +
geom_errorbar(aes(ymin = cost_agb - SE_agb, ymax = cost_agb + SE_agb),
width = .2, position = position_dodge(0.9)) +
scale_fill_manual(
labels = c('5 years', '7 years'),
values = c("#CAF0F8", "#00B4D8")) +
geom_text(data = letras, aes(x = tratamento, y = y, label = label), size = 3.5, color="black", vjust = -0.5) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 240), breaks = c(0,50,100, 150, 200))+
labs(fill = "",
y = expression("CE AGB (" * "\u0024" * ".Mg"^{-1} * ")"),
x = "") +
theme_classic()+
theme(
axis.line.x = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
axis.line.y = element_line(color = "black", linewidth = 0.8, linetype = "solid"),
axis.ticks = element_line(color = "black", linewidth = 0.5),
axis.text.x = element_text(color="black", hjust = .5, size=10),
axis.text.y = element_text(color="black", size=10),
legend.position = 'none',
axis.title.y = element_text(color="black", size=10))
fig_cost2
show code
fig_cost <- cowplot::plot_grid(fig_cost1, fig_cost2, ncol = 1,
align = "v", axis = "l",
labels = c("(a)", "(b)"), label_size = 8,
label_x = 0.05, label_y = c(1,1.1))
fig_cost
show code
# ggsave("fig_cost.png", plot=fig_cost, width=90, height = 120 ,unit="mm")(iii) Table Cost-effectiveness
show code
fig_data_cost_raw |>
dplyr::mutate(dplyr::across(where(is.numeric), ~ round(.x, 2))) # A tibble: 10 × 6
ano tratamento cost_cc SE_cc cost_agb SE_agb
<dbl> <ord> <dbl> <dbl> <dbl> <dbl>
1 5 "Nuclei\n25%" 23.3 9.43 133. 69.9
2 5 "Nuclei\n50%" 15.5 1.08 46.4 4.84
3 5 "Planting\n100%" 7.07 1.09 25.6 3.28
4 5 "Strips\n25%" 12.8 1.8 44.6 5.95
5 5 "Strips\n50%" 14.5 3.44 43.3 8.74
6 7 "Nuclei\n25%" 10.5 4.82 58.6 26.0
7 7 "Nuclei\n50%" 6.26 0.94 24.4 2.88
8 7 "Planting\n100%" 5.35 0.11 16.6 1.47
9 7 "Strips\n25%" 7.52 1.37 22.6 2.66
10 7 "Strips\n50%" 6.63 1.92 26.3 7.65